Emergent smectic order in simple active particle models

Novel"smectic-P"behavior, in which self-propelled particles form rows and move on average along them, occurs generically within the orientationally-ordered phase of simple models that we simulate. Both apolar (head-tail symmetric) and polar (head-tail asymmetric) models with aligning and repulsive interactions exhibit slow algebraic decay of smectic order with system size up to some finite length scale, after which faster decay occurs. In the apolar case, this scale is that of an undulation instability of the rows. In the polar case, this instability is absent, but traveling fluctuations disrupt the rows in large systems and motion and smectic order may spontaneously globally rotate. These observations agree with a new hydrodynamic theory which we present here. Variants of our models also exhibit active smectic"A"and"C"order, with motion orthogonal and oblique to the layers respectively.


Introduction
Active matter, by which we mean out of equilibrium systems that locally convert energy into directed motion, can potentially exhibit new phases not present in equilibrium systems. Even active phases that do have equilibrium counterparts often behave quite differently from them, in both their fluctuations [1,2,3,4,5,6,7] and flow properties [8,9,10,11,12,13].
All of the work just described and cited dealt with phases of active particles with orientational order, but no translational order [14,15]. More recently, attention has turned to translationally ordered phases [16,17,19,18,20,21,22]. One class of such phases is "active smectics" [23,24], which are out of equilibrium analogs of smectic liquid crystal phases. In these phases, in addition to orientationally ordering, the particles spontaneously form regularly spaced, liquid-like layers. That is, translational symmetry is broken in only one direction, while the systems themselves are two or three dimensional.
This theoretical work focused on active analogs of the smectic-A phase, for which particle alignment is perpendicular to the layers. These theories [23,24] predicted in particular that smectic order is long-ranged in d = 3 dimensions and quasi-long ranged (i.e., power-law correlated) in d = 2, in strong contrast to the equilibrium case, in which smectic order is only quasi-long ranged [26] in d = 3 and completely destroyed by thermal fluctuations (i.e., short-range correlated) in d = 2 [25]. This phenomenon -that is, activity stabilizing order in low spatial dimensions in which it is absent in equilibrium-is similar to the "violation" of the Mermin-Wagner theorem that occurs in polar active fluids in two dimensions [1,2,3].
Prior to the work we report here, none of the predictions made in [23,24] had been confirmed either in experiments or simulations. For active particle systems, few striped patterns have been reported so far [27,28,17,29], and none have been studied quantitatively as smectic phases. (See, however, the shaken rice grain experiments of Narayanan et al. [30].) In this paper we report the generic emergence of active smectic order in simulations of very simple self-propelled particle models, in which the alignment interactions of the familiar Vicsek [31] model are supplemented by short-ranged repulsive interactions between particles. Such models have been studied [7,32,22], but not in the high density regimes we consider here. Simulating both this modified Vicsek model and an apolar (i.e., head-tail symmetric) version of it in d = 2 space dimensions [6], we observe the development of quasi-long-ranged smectic order over a range of length scales in part of the orientationally-ordered phases exhibited by these models (cf. Fig. 1). When the repulsion is anisotropic (i.e., dependent on the angle between the direction of particle motion and the inter particle relative position vector), we find also smectic-A (in which the layers are orthogonal to the particle alignment) and smectic-C order, in which the layers make an acute angle with the particle alignment ( Fig. 2 a,b,c). Both of these configurations also occur in equilibrium. But the active smectic configuration we most frequently observe, particularly with isotropic repulsion, is a new state of active matter with no equilibrium counterpart that we call"smectic P", in which the particle alignment is primarily parallel to the layers (Fig. 1b,c). A pictorial summary of the polar and apolar smectic A and P phases is given in Fig. 3.
We deliberately refrain from calling the regions of smectic order we have found phases, since in all cases studied, quasi-long-ranged smectic order disappears beyond some model-dependent length scale. This disappearance of order is explained by comprehensive phenomenological theories of polar and apolar active smectics that we present here. These theories predict that symmetry-allowed nonlinearities in the polar case destroy smectic order at the longest length scales, and that a long but finite wavelength "undulation" instability likewise disorders the apolar case. Both predictions agree with our simulations.
Finally, we also report a number of spectacular collective dynamical phenomena, such as spontaneous global rotation of both orientational and smectic order for some polar smectic P, and large-scale swirling structures in polar smectic A, which await theoretical elucidation.

Simulation models
We consider two classes of models: In the polar model, particles move along an intrinsic heading vector, which they align ferromagnetically with that of their neighbors. In the apolar model, particles align their axes with those of their neighbors, and move with equal probability in either direction parallel to this axis. In both classes, we add to this alignment interaction (which is in competition with noise) pairwise repulsion between neighbors. We consider both isotropic repulsion, in which the repulsion does not depend on the angular position of neighbors with respect to the particle's intrinsic axis, and anisotropic repulsion, in which it does.
More specifically, in the spirit of the Vicsek model [31], we consider point particles moving at constant speed v 0 . At discrete timesteps, the position r i of particle i is moved along the unit vector u i (t + 1) = (cos θ i (t + 1), sin θ i (t + 1)) T : where χ i (t) ∈ [− π 2 , π 2 ] is an angular white noise drawn from a uniform distribution, σ is a parameter setting the strength of the angular noise, and A i and R i are respectively "alignment" and "repulsion" vectors obtained from averages over neighbors defined as those particles located within a distance r int ., with the default interaction range set to unity (i.e. r int = 1, see Appendix F). The same unit range was chosen for both interactions for simplicity. We have checked that taking different ranges does not qualitatively change the model's behavior.
The two models are distinguished by their alignment interaction: where the N i neighbors j include particle i, (t) = 1 for the polar model, while for the nematic model (t) is randomly chosen to be ±1 at each step with equal 50% probability. The repulsion vector R i is identical for both cases: it is the average over all neighbors of a pairwise force alongr ji , the unit vector pointing from particle j to i. The magnitude of this repulsive force is the same for all neighbors in the isotropic case, while, in the anisotropic case, it depends on φ ji = acos(u i ·r ji ), the angular position of particle j relative to the axis of particle i Here we only consider γ = 0 and γ = π 2 . For γ = 0, repulsion is stronger ahead of and behind the particles, whereas for γ = π 2 it is stronger towards the left and right.

Regions of smectic order
For β = 0, the polar and apolar models defined above respectively reduce exactly to the Vicsek model [31] and the active nematics model of [6], which both have two main parameters: the noise strength σ and the global mean number density ρ 0 . Thus, β, which controls the repulsion strength, is the only new parameter. We have explored systematically the three-parameter space of our models. As in the repulsion-free case, an orientational order-disorder transition is found upon increasing σ and/or decreasing ρ 0 , in both the polar and the apolar case. Note that the coexistence phase of this transition, best described in the liquid/gas framework [33], is essentially unobservable at such high-densities. Significant smectic order is only found for relatively high densities (ρ 0 8), inside the orientationally-ordered phase (Fig. 1a). We quantify this smectic order with smectic order parameters S n , which are defined as S n ≡ I(q n , t) , (2.5) where I(q, t) ≡ |ρ(q, t)| 2 /N 2 , with N = ρ 0 L 2 the total number of particles, and q n ≡ nq 0ẑ , with q 0 ≡ 2π a l . Here ρ(q, t) = d 2 r ρ(r, t)e −iq·r .
is the spatial Fourier transform of the density, n is a non-zero integer,ẑ is the mean normal to the layers, and a l the layer spacing, which is typically of the order of the  (red oblique crosses) falling perfectly on top of S 1 for L L c , as predicted by the hydrodynamic theory. Above L c scaling breaks down due to the undulation instability. (g) The polar model with ρ 0 = 10, β = 0.05, σ = 0.007 shows similar weak algebraic decay at small L L N L ≈ 60. However for L L NL the results for S 1 , S 1/4 2 , and S 1/9 3 start to diverge, and S 1 decays faster than algebraically. interaction range, see Appendix E. We hereafter restrict ourselves to high average densities (ρ 0 ≥ 8) and low speeds (v 0 = 0.2 − 0.3), and vary σ and β. We use square domains of linear size L with periodic boundary conditions, and consider only "zerowinding number" smectic configurations (see Appendix G for details).
The polar and apolar models, whether with isotropic or anisotropic repulsion, share qualitatively similar phase diagrams in the (β,σ) plane, whose main features are shown in Fig. 1d for the apolar model with isotropic repulsion. When both repulsion and noise strength σ are sufficiently low, the system develops orientational order. Inside this region, a line starting at β = σ = 0 and finishing at β max close to the orientational orderdisorder transition delimits a domain in which smectic order appears over a large, but finite, range of length scales (cf. Appendix E, Fig. E1). The location of the boundaries of these regimes do not change much with system size L, but the maximal levels of global smectic order generally decrease with increasing L (whereas local smectic order remains strong).
The type of smectic order observed depends on the model (Table 1). For isotropic and anisotropic repulsion with γ = π 2 , a novel type of smectic order emerges that is not observed in equilibrium [34]. We call this new type of smectic order, in which the particle axes are parallel to the layers (Figs. 1b,c & 3b,d), "smectic P".
Since repulsion between particles favors translational order, anisotropic repulsion favors anisotropic translational order (i.e., layering). This implies, in particular, that when repulsion is stronger perpendicular to the direction of particle motion (as is the case for case γ = π/2), the smectic P configuration just described is preferred.
In our case, isotropic repulsion leads to the smectic P phase, which also occurs, for the reasons just given, in the explicitly anisotropic case γ = π/2, for which repulsion is stronger to the left and right of the direction of motion of the particles. This is because, in the isotropic case, particles do not move away from those ahead of or behind them as effectively, since this requires speeding up or slowing down, which in the simplest model with fixed speed they cannot do. We have confirmed that the P phase is also favored for active particles models with variable speed [35], in which the speeds of the particles are allowed to fluctuate around the preferred speed v 0 due to the action of the various forces. In general, the additional anisotropy linked to the propulsion direction of self-propelled particles effectively mimics the effect of weaker repulsion in those directions.
For γ = 0, repulsion is strongest in the front and back of particles, and one might naively expect smectic-A order. However, in the apolar case, this is in fact never observed. Instead, the left/right symmetry is spontaneously broken, and the particles' axes tip away from the layer's normal: a smectic C with |Φ − Φ S | 0.2π, where Φ and Φ S denote respectively the angles of the particle axes and the layer normal (cf. Fig 2  d). Only in the polar model, after a possibly long transient during which patches of smectic-C order compete and form locally chevron-like structures, Fig. 2b), does the system eventually reach a smectic-A state (Fig. 2c). Table 1. Type of smectic phase (A, C, P) and large-scale phenomena observed for polar and apolar models with isotropic and anisotropic (γ = 0, π 2 ) repulsion.

Quasi-long-range order and large-scale phenomena
In all cases, defect-less smectic configurations are only observed in sufficiently small (albeit still quite large) systems (see Supp. Mat. Movies 1 & 2). Even in simulations that start from carefully prepared "perfect" initial configurations, spontaneous nucleation of dislocations and/or large-scale undulations occur beyond some model-dependent lengthscale, which we call L c (L NL ) for the apolar (polar) case; the meaning of these names is explained in the "Hydrodynamic theory" section below. For L less than this length scale, S n shows either no decay, or very slow algebraic decay with L; the latter is what we mean by quasi-long-range smectic order. For much larger systems, S n decreases much faster with L ( Fig. 1f,g). For L L c in the apolar case, the system shows undulations of the smectic-P layers (Fig. 1e). The wavelength of the undulations is well defined asymptotically and typically large (∼ 100 − 200). Note that this wavelength is of order of L c . This undulation pattern is not quite steady: its "arches" slowly move as indicated by the arrows in Fig. 1e and Supp. Mat. Movie 3. This accompanied by the constant nucleation of dislocations in the sheared regions between arches. For L L NL in the polar case, no undulation instability is observed. At large scales, polar smectics-P layers show traveling fluctuations on all scales. In contrast to the apolar case, these fluctuations never cohere into a regular moving pattern that is stationary in a co-moving frame. Instead, the deviations continue evolving randomly in time, and do not develop a single typical wavelength. Note that our theory predicts that fluctuating modes with the largest wavelength are expected to have the largest amplitude, which may produce the impression of a transient regular pattern with wavelength of the order L/2 (Supp. Mat. Movie 2).
We finally report on two spectacular phenomena observed with anisotropic repulsion. For polar smectics P with γ = π 2 , we observe a spontaneous breaking of the left/right symmetry in the form of global rotation of both the particle axes and the smectic layers. For small β values, at fixed system size, rotation may be intermittent and may change sign in time ( Fig. 2d and Supp. Mat. Movie 4). But for sufficiently large systems, even an initially prepared defect-free smectic-P configuration starts rotating either clockwise or counterclockwise with a well-defined, steady angular velocity ω which  increases strongly with β (Fig. 2e). Our observations of the intermittent, slow-rotation regimes reveal that the polar orientation angle Φ is 'driving' rotation: the angle of particle axes changes first, with Φ S , the angle of the layer's normal, fixed, lagging behind. The stress thus induced on the smectic layers then generates dislocations, Φ S slips, and the defects annihilate. Since Φ continues rotating, this sequence continuously repeats itself. We note that we do not observe hopping of individual particles between the layers; instead the observed global rotation corresponds to breaking of layers and reattaching to neighboring layers, which leads to formation of dislocation lines travelling through the system ("defect waves", Supp. Mat. Movie 4). For the polar smectic-A, another spectacular phenomenon sometimes occurs for the larger β values within the smectic region: large-scale rotating "swirls", made of a spiral, flower-like arrangement of smectic layers (Fig. 2g, and Supp. Mat. Movie 5). In a periodic domain, a number of clockwise and counterclockwise swirls may emerge, which experience some effective repulsion, leading to stable configurations of equal number of '+' and '−' swirls that coexist with the swirl-free one.  These systems exhibit smectic-P order for system sizes L L c , and instability for L L c . The critical system size for the parameters used here is L c ≈ 200. In contrast to both polar and apolar active fluids with orientational, but no translational, order, we observe no "giant number fluctuations": that is, we do not observe ∆n 2 ∝ ∆n 2 α with α > 1; rather, we find conventional behavior (i.e., α = 1). Parameters: ρ = 8, v = 0.3, σ = 0.01, β = 0.1

Hydrodynamic theory
We now demonstrate that most of the above results are well accounted for by hydrodynamic theories of active smectics. The variables of the hydrodynamic theory are the field u(r, t) giving the displacement of the layers perpendicular to their unperturbed, uniformly spaced, parallel locations, and the number density ρ(r, t). The latter is a slow hydrodynamic variable because of number conservation, while the former is the Goldstone mode associated with the spontaneous breaking of continuous translational symmetry in the direction perpendicular to the layers. The theory consists of a closed set of stochastic partial differential equations for the time evolution of u(r, t) and ρ(r, t), containing all terms at lowest order in a gradient expansion that are consistent with the symmetries of the underlying dynamics and the broken symmetry state. Since we have two different broken symmetries (polar and apolar smectic P), we have two sets of hydrodynamic equations.

Apolar smectic P
In two dimensions, there is no symmetry difference between the smectic-P order found here and the apolar active smectic-A phase treated in [23]: both are symmetric under separate inversions about the z axis (the normal to the smectic layers) and about the x axis (along the layers). Hence the equations for this case are the same as those developed in [23] for apolar active smectic A. They read where δρ ≡ ρ − ρ 0 , and f u and f ρ are Gaussian, zero-mean, white noises with variances We set the cross-correlation f ρ (r, t)f u (r , t ) = 0; corrections to this can be shown to be "irrelevant" in the renormalization group sense of having no effect on the large distance, long time behavior of the system. Here B, D ux , K, D ρx , D ρz , C, C x , C z , ∆ u , ∆ ρx , and ∆ ρz are all phenomenological parameters that cannot be determined by symmetry arguments, but must be deduced from experiments, simulations or a detailed kinetic theory.
In an equilibrium smectic, the D ux term in (3.1) is forbidden by rotation-invariance of the free energy. It is, however, permitted here [36] simply because rotation-invariance at the level of the EOM, which is all one can demand in an active system, does not rule them out. The physical content of this, term is that layer curvature produces a local vectorial asymmetry which must lead to directed motion of the layers, as this is a driven system. A similar term occurs in single membranes with "pumps" [36,37].
The B term simply acts to keep the smectic layers equidistant, while the K term is a bend modulus that tends to keep them straight. The rest of the terms in (3.1) and (3.2) are present in equilibrium, although the equilibrium requirement [38] that Cx Cz = Dρx Dρz does not hold in our non-equilibrium system. There is no a priori symmetry argument that determines the sign of D ux . As shown in [23] when D ux > 0, smectic layer fluctuations are suppressed, leading to quasilong-range order in d = 2, in contrast to the short-ranged order found in equilibrium [25]. When D ux < 0, the perfect smectic state with constant u(r, t) and ρ(r, t) = ρ 0 is unstable against "undulations", in which the layers wiggle in unison. This is easily seen by Fourier transforming the (linear) equations of motion in space. For wavevectors q with q z = 0, ρ and u decouple, leading to the growth rate ν(q We show in Appendix A.1 that allowing q z = 0 does not change this result: the instability always grows fastest along x. This is precisely what our simulations show: undulations appear and grow along the layers, but only in systems whose extension L x along this axis is large enough. This minimal size is of the order of the scale of the observed undulations, which we identify with: Note that the wavelength of the undulations is of order L c , but not necessarily equal to it, as it will ultimately be determined by nonlinear, saturating terms neglected here. We will also use Eq. (3.5) to define a "critical length scale" in the stable case D ux > 0, as the scale beyond which active tension effects (which then suppress fluctuations) become important.
(2) L NL Figure 5. Regions distinguished by the scaling behavior of smectic order parameters S n with system size L and active tension D ux . Red lines indicate the lengthscale L c defined in the main text. In region I, defined by L < L c ≡ 2π K/|D ux | and L ξ x , in the apolar case, S n ∝ L −n 2 η I ; in the polar case, S n 's do not decay. In region II, defined by L < L c and ξ x L ξ z , both the apolar and polar smectic P obeys S n ∝ L −1−η I n 2 n −2 . In region III (L < L c and L ξ z ), both the apolar and polar smectic P obey S n ∝ L −2 n −6 : purely short-ranged smectic correlations, as explained in the main text. In region IV (L L c and D ux > 0) S n ∝ L −n 2 η IV for both the apolar and polar cases. In region V, which is L L c and D ux < 0, both the apolar and polar smectics P are unstable, and we expect S n ∝ L −2 . Region VI, in which L L NL , only exists for the polar case. Here nonlinear effects become important, and may induce dislocations that destroy the smectic order. If they do, then we expect S n ∝ L −2 in this region as well. The horizontal dashed lines above and below the L-axis are the path followed by increasing L at fixed parameters in our simulations of respectively the apolar model (1, data shown in Fig. 1f) and the polar model (2, data shown in Fig. 1g).
We now analyze noise-induced fluctuations in our theory. We focus on the case in which L c is large, the positional noise ∆ u is small, and the density noises ∆ ρ(x,z) are large, all of which, as we argue in Appendix A.2, are satisfied in our simulations. Here we simply sketch the derivation of, and summarize, our results; details are given in Appendices A. 3

and A.4.
We begin by spatiotemporally Fourier transforming the equations of motion (3.1) and (3.2) and solve the resultant linear algebraic equations for the transformed fields u(q, ω) and ρ(q, ω) in terms of the random forces f u (q, ω) and f ρ (q, ω). Autocorrelating, and using (3.3) and (3.4) gives an expression for C uu (q, ω) ≡ |u(q, ω)| 2 . Integrating this over all frequencies ω gives the equal time correlation function C ET uu (q) ≡ |u(q, t)| 2 , from which we can in turn calculate equal-time real space correlation functions.
The correlations of u we thereby obtain are related to the smectic order parameters. Recall that S n ≡ I(q n , t), where I(q, t) ≡ |ρ(q, t)| 2 /N 2 , with N = ρ 0 L 2 the total number of particles, and q n ≡ nq 0ẑ . These S n are sensitive to fluctuations of the displacement field u since a translation by ∆u clearly changes the complex phase of ρ(q n , t) by nq 0 ∆u. We show in Appendix A.4 that for the apolar model: where |∆u(r)| 2 ≡ |u(r + r , t) − u(r , t)| 2 can be easily obtained by Fourier transforming C ET uu (q) back to real space, and the w n 's are uninteresting constants given in Eq. (A.70). This calculation, given in detail in Appendix A.3, reveals the existence of two other important lengths in addition to L c : The resultant behavior of S n with system size exhibits numerous crossovers, which delimit five different regions of the (L, D ux ) plane, as summarized in Fig. 5: Here the order parameters S n decay algebraically with L: S n ∝ L −n 2 η I , with η I a non-universal exponent that is a monotonically increasing function of the noise strengths ∆ ρx and ∆ ρz .
Region II, in which L L c and ξ x L ξ z , is characterized by rather fast algebraic decay of S n according to S n ∝ L −1−η I n 2 n −2 .
Region III, defined by L L c and L ξ z , shows even faster decay with L: S n ∝ L −2 n −6 . This corresponds to purely short-ranged smectic correlations. the system breaks up into decorrelated smectic regions of fixed size whose number scales as L 2 ; their contributions to the global order parameters S n simply add randomly, leading to the scaling S n ∝ 1 Nr ∝ 1 L 2 . Region IV, in which L L c and D ux > 0, exhibits quasi long-ranged smectic order with S n ∝ L −n 2 η IV , where η IV is very close to, but different from η I (see Eq. A.74) Region V is defined by L L c and D ux < 0. Here, the system is unstable against undulations, and we expect S n ∝ L −2 .
Note that varying L, while keeping parameters fixed, will not necessarily explore all of these regions. For example, a system with a sufficiently large in magnitude, negative D ux (specifically, large enough that L c ξ x , which would not have to be very large if the positional noise ∆ u is small, since then ξ x is large) will follow the horizontal locus labeled "(1)" in Fig. 5, and will pass directly from regime I to the unstable regime V. This is precisely the scenario we observe in our simulations of the apolar model. Our numerical results show that for L < L c , S 1 (L), S 2 (L) 1/4 , and S 3 (L) 1/9 coincide (Fig 1f), which is in full agreement with our theoretical prediction in region I, in which S n satisfy (3.8) For L L c , S 1 (L), S 2 (L) 1/4 , and S 3 (L) 1/9 depart from each other (Fig 1f), and the undulation instability appears (Fig. 1e), which implies the system is in the unstable region V. As mentioned earlier, the slow drift "arches" of the undulations in our simulations generates dislocations (Fig. 1e), which destroys the smectic order. Therefore we expect to see a much faster algebraic decay of S n with L in region V than in region I. This is, indeed, what we see, as shown in Fig. 1f. Note also that the weak algebraic decay of S n for L L c observed in our simulations is not that predicted in [23]. That prediction applied in region IV in Fig. 5, which is the region of longest length scales in the stable case (D ux > 0).

Polar smectic P
Polar smectics P do have a different symmetry than polar smectic A, because the mean motion parallel to the layers (i.e., along the x-direction) breaks the x → −x symmetry. Adding to the original apolar model (3.1) and (3.2) all x-inversion symmetry-breaking terms that are "relevant", the hydrodynamic equations read ‡: The terms in these equations that are identical to those in (3.1) and (3.2) have the same physical origin as they did there. The new terms (e.g., the g term), are present here, but not in the apolar case, because they violate x → −x symmetry, which is not present in the polar case, but is present in the apolar case. The noises f u and f ρ have the same statistics (3.3) and (3.4) as before. As in the apolar case, the trivial ordered solution undergoes a long-wavelength instability at zero noise when D ux < 0. The dispersion relation for q z = 0 modes now reads x . Hence the instability (which again is strongest for q alongx) has the same spatial structure as before, but propagates dispersively due to the g term, with phase velocity v p (q x , q z = 0) = gq 2 x . We thus expect (assuming as usual that nonlinear saturation does not modify the mode structure) the most unstable mode to propagate at a speed v c given by . For D ux > 0, modes of any wavevector q are now stable, but each propagates along the layers with the phase velocity v p (q x , q z = 0) = gq 2 x . We observe such propagation of fluctuations in our simulations (Supp. Mat. Movie 2).
Linearizing Eqs. (3.9) and (3.10), we can calculate the behavior of the smectic order parameters S n with L exactly as we did in the apolar case. Details are given in Appendix B; the result is that the active polar smectic P, in this linear approximation, ultimately has the same crossover structure as that illustrated for the apolar case in Fig.  5. The only difference is that η I = 0, which implies that in region I, the S n 's, rather than falling off algebraically with L, are essentially constant.
The above discussion was based entirely on the linear approximation. As mentioned earlier, the nonlinearities in Eqs. (3.9), (3.10) are relevant (technically, "marginal") in d = 2. We therefore expect that there exists a nonlinear length scale L NL , growing very rapidly (specifically, like exp(A/∆ ν ), with A a non-universal constant and ν a universal exponent that we have not yet determined) with decreasing noise strength ∆, where by ∆ we mean a suitably weighted average of the noise strengths ∆ u,ρx,ρz , beyond which our theory is not valid. What happens beyond L NL remains an open question, which can only be answered by a full renormalization group analysis. Since such an analysis is quite formidable, we restrict ourselves to plausible speculation. The relevant nonlinearities may make the system much softer (i.e., fluctuate much more strongly) at long wavelength, as in other similar dynamical problems (e.g., the two dimensional KPZ equation [39]). Then this softness will almost certainly lead to the unbinding of dislocations at the longest length scales L L NL , even when D ux is positive. This would imply that all polar smectic-P's would be disordered at the longest length scales. Therefore, we expect both a breakdown of scaling law (3.8) and a much faster algebraic decay of S n (L) with L for L L NL . This adds another region to Fig. 5: Region VI : for polar smectic P with L L NL , nonlinear effects become important, which may induce dislocations and destroy smectic order. We expect S n (L) ∝ L −2 .
This new region (VI) may not exist for all polar smectics P. In active polar smectics A [24], nonlinear effects of the same type (i.e., marginal in d = 2) as those found here destroy smectic order at long length scales in some, but not all, regions of parameter space. Whether this can happen for polar smectics P remains an open question.
These predictions and speculations are in good agreement with our simulations, as illustrated in Fig. 1g, in which we again plot S 1 (L), S 2 (L) 1/4 , and S 3 (L) 1/9 . They lie on top of each other for small system size, as predicted above for region I and IV. For large L they depart from each other and decay much faster: we interpret the system size at which this happens as L NL . That is, these simulations are following the locus labeled "(2)" in Fig. 5. Apparently, the active tension D ux is so large that the locus does not enter into regions II and III. However, unlike in the apolar case, D ux appears to be be positive here, since we observe no undulation instability. Although, we emphasize, this is not a universal property of polar smectic-P order: it is entirely possible that microscopic models different from ours may exhibit an undulation instability.

Summary and discussion
We have shown that active particle models which have both repulsion and alignment generically exhibit smectic configurations for sufficiently large densities. This occurs in both apolar models, in which the particles are equally likely to move in either direction parallel to their body axes, and in polar models, in which there is a preferred sign of motion along the axis.
These results were obtained with Vicsek-style models, but we have checked that our conclusions also hold for continuous-time Langevin models of the type studied in [35]; ergo, they are not artifacts of the discrete-time updating. Depending on the symmetry of the particles and alignment, and on the anisotropy of the repulsion, smectic order and various dynamical large-scale phenomena emerge, as summarized in Table 1.
We have found that a new type of smectic arrangement is most common in our models: a state we call "Smectic P", in which particle alignment (and motion) is primarily parallel to the layers (Fig. 1b,c). This smectic-P order is unique to active particles: states with the particle axes along the smectic layers have never been seen in equilibrium.
In both the apolar and the polar cases, we observe weak algebraic decay of smectic order parameter S n (L) with system size L for L smaller than some crossover length, and much faster algebraic decay for larger L. In the apolar case we also observe an undulation instability in large systems, while in the polar case we do not.
To understand these phenomena, we have modified the hydrodynamic theories of active smectics A introduced in [23,24] to treat both apolar and polar smectics P. The theory of apolar smectics P proves to be the same as that for apolar smectics A; indeed, we interpret the instability we see as precisely the instability predicted in the "negative active tension" case in [23]. We have, however, extended this theory of the apolar case to treat the regime of system sizes L smaller than the instability length L c . We find that two additional important length scales smaller than L c can also appear for small values of the "active tension" D ux , which is the most important fundamentally non-equilibrium parameter in the hydrodynamic theory. This theory leads to the prediction of the five distinct regimes of behavior in (L, D ux ) plane illustrated in Fig. 5.
The behavior that we observe in our simulations of the apolar case, as just summarized above, is entirely consistent with this hydrodynamic theory for a system with sufficiently negative active tension D ux , such that, with increasing system size L, our system moves along the locus labeled (1) in Fig. 5. Note that this instability is not inherent to all apolar systems. Some models that have the same symmetries as ours, but that differ in detail, may exhibit stable apolar phases. The hydrodynamic theory predicts that such a stable apolar system will exhibit algebraically slow decay of smectic order out to arbitrarily large system size, at least for sufficiently small noise.
The hydrodynamic theory of the polar phase predicts rather different behavior. In addition to the five regions of the (L, D ux ) plane, a sixth region appears, which includes all systems larger than yet another characteristic length L NL (Fig. 5). We believe dislocations will always appear in this region, and make smectic order short ranged. Our numerical results are consistent with this prediction for a system with large positive active tension D ux .
Note that, just as not all apolar systems need be unstable, nor do all polar systems need to be stable. Some polar models, different from ours, will probably exhibit an undulation instability of the type we observe in the apolar case, although, in the polar case, this instability be propagative.
We did not observe "giant" number fluctuations, which are a signature of orientationally ordered but translationally disordered phases [1, 4, 7, 30] (Fig 4). But these simulations, because of their intrinsic difficulty, were performed at sizes never much larger than the characteristic lengths L c or L NL . We cannot exclude, and indeed believe, that giant number fluctuations, as well as other exotica generically present in orientationally-ordered flocks, will exist in very large active systems showing local smectic order and long-ranged orientational order.
Despite the general success of the hydrodynamic theories developed here in explaining many of the phenomena we observe in our simulations, some of them remain mysterious. In particular, the spontaneous rotation of the smectic layers (Fig. 2a,b, Supp. Mat. Movies 4 & 5) is a type of spontaneous chiral symmetry breaking that is beyond the scope of the theories presented here, which assume an achiral steady state. We hope to develop a hydrodynamic theory of such symmetry breaking in future work.

Acknowledgments
We thank F. Ginelli and S. Ramaswamy for lively and illuminating discussions which took place within the Advanced Study Group "Statistical physics of collective motion", generously supported by the Max Planck Institute for the Physics of Complex Systems in Dresden. J.T. also thanks the Aspen Center for Physics, Aspen, Colorado; the Isaac Newton Institute, Cambridge, U.K., and the Kavli Institute for Theoretical Physics, Santa Barbara, California; for their hospitality while this work was underway. He also thanks the US NSF for support by awards # EF-1137815 and 1006171; and the Simons Foundation for support by award #225579. L.C. acknowledges support by the National Science Foundation of China (under Grant No. 11474354) Appendix A. Hydrodynamic theory predictions for the apolar active smectic P phase Appendix A. 1

. Eigenfrequencies and Instability threshold
We begin by Fourier transforming the equations of motion (3.1), (3.2) in the main text: where we have defined: The two eigenfrequencies ω ± of these equations of motion are, as usual, those values of ω that allow them to have non-zero solutions for δρ and u when the forces f u and f ρ on the right hand sides are set to zero; i.e., they are the solutions of the eigenvalue equation: which are: Stability requires that the imaginary part of these eigenfrequencies to be negative. The real part is readily seen to be zero. This will be the case if and only if two conditions are satisfied: Squaring the first of these conditions (A.8) implies which can be reorganized to read: Using our definitions (A.3), (A.4), and (A.5) for Γ u (q), Γ ρ (q), and Γ ρu (q) in this expression gives, after gathering terms, where we have defined A ≡ BD ρz − CC z and F ≡ BD ρx − CC x + D ux D ρz . Note that A and F will both be > 0 for sufficiently small C, provided that D ux is not too large and negative. In this case the left hand side of (A.12) is a monotonically increasing function of q 2 z , and, hence, has its minimum at q z = 0. Hence, it is at q z = 0 that the condition (A.12) will first be violated. Furthermore, this violation will first occur with decreasing q x when q x = −D ux /K = q max , where q max is the critical value of q x defined in the main text. (Recall that here we are discussing the case D ux < 0.) Thus, our claim in the main text that the first instability with increasing system size occurs at q z = 0 , q x = q c follows provided that our second condition (A.9) is not violated at a larger q. But noting that Γ u (q) and Γ ρ (q) are both > 0 for any q x > q c (since then D ux q 2 x + Kq 4 x > 0, and we also know that B and both D ρx and D ρz are > 0), it is clear that (A.9) can not be violated before (A.8). This completes our demonstration that the first instability with increasing system size occurs, in the apolar case, at q z = 0 , q x = q c .

Appendix A.2. Fourier space u − u correlation functions
The linear algebraic equations (A.1) and (A.2) are easily solved for u(q, ω) and δρ(q, ω). The result for u(q, ω) is where ω ± are the two eigenfrequencies found above. From equation (A.6), we can read off and both of which will prove useful later. We also note that ω ± are purely imaginary, which of course implies that ω * ± = −ω ± . Using this last fact, and our expression (A.13) for u(q, ω), we immediately get for the autocorelations which is obviously born to be integrated by complex contour techniques. Before doing so, however, we need the spatio-temporally Fourier transformed autocorrelations of the noises |f u (q, ω)| 2 and |f ρ (q, ω)| 2 . These are easily read off from the real space correlation functions (3.3),(3.4) in the main text: Using these in (A. 16) gives Integrating (A.19) over all frequencies ω by obvious complex contour techniques then gives, after a little algebra, the equal time correlation function: which can be simplified using our earlier expressions (A.14) and (A.15) for the sums and products of the eigenvalues, yielding, after a slight rearrangement of terms, our final result: where we have defined , (A. 23) and . (A.24) The first term C ET (1) uu in (A.21) exhibits very different behavior for wavevectors with q x q c ≡ |D ux |/2K and q x q c . Note that for D ux < 0, q c is the wavevector of maximum instability described in the main text. For D ux > 0, q c is simply a crossover wavevector between the two regimes that we'll now describe. For q x q c , we can, by the definition of q c , neglect the D ux term in Γ u . Doing so, we immediately see that the first (i.e., the Γ u ) term in the correlation function (A.20) has the same form, for wavevectors q with q x q z , as the full, equilibrium [25] C uu (q). It is given by where we have defined the "smectic penetration depth" λ ≡ K/B and B ≡ B − CC x /D ρx . Note that this is much larger than q −2 for q x q z ; hence, this range of wavevectors dominates the contribution of C ET uu (q) to the real space fluctuations for system sizes L L c . Note also that C ET uu (q) is, for small D ux , extremely anisotropic: for q x q c , it scales like q −2 for q z λq 2 x , and like q −4 for q z λq 2 x . For q x q c , it is ∼ ∆u Bq 2 for q z D ux /Bq x , while it is ∼ ∆u Duxq 2 for q z D ux /Bq x . In short, it is always much larger, for a given |q|, for small q z ; i.e., for directions of wavevector q near the plane of the smectic layers. It also changes its behavior dramatically between the range q x q c , in which it scales like q −2 for all directions of q, and q x q c , where it scales like q −2 for q z λq 2 x , and like q −4 for q z λq 2 x . The second term in (A.21) has neither this anisotropy, nor this sensitivity to whether q x q c or q x q c . To see this, note that the ratio Cq 2 z Γ u (q) is bounded above by C B , while the sum Γ u (q) + Γ ρ (q) scales like q 2 for all directions of q, even when D ux is small. (Note that this is true even though we are considering q x q c , which is, by definition, the regime of wavevectors in which the Kq 4 x term dominates the D ux q 2 x term, because the D ρx q 2 x term in Γ ρ still dominates the Kq 4 x term, since D ρx , unlike D ux , need not be small when the activity is small.) Indeed, this second term scales like q −2 for all directions of q. We'll see in a moment that, as a result, this term leads to essentially isotropic algebraic decay of smectic order parameter correlations at all length scales.
Note also that only this second term grows with increasing density noise ∆ ρ(x,z) (since the first term depends only on the positional noise ∆ u ).
The above results are completely general. We will in what follows frequently consider the limit of large L c , small positional noise ∆ u , and large density noises ∆ ρ(x,z) .
All three of these criteria are satisfied in our simulations: empirically, L c ∼ 100a; the positional noise ∆ u should be proportional to the mean squared velocity fluctuations perpendicular to the layers, which, since the motion is primarily parallel to the layers, should scale like our angular noise strength parameter σ 2 ∼ 10 −4 . Finally, the density noises ∆ ρ(x,z) are large, because our algorithm introduces an O(1) noise in the step that randomly, with equal probability, makes the particle move either forward or backward along the direction selected. This will not, in the small or zero σ limit, contribute to the displacement noise ∆ u , since in that limit all of the motion is along the layers. But random statistical fluctuations in the number of particles moving left or right within a layer will clearly lead to fluctuations in the density; that is, to appreciable ∆ ρ(x,z) . This combination of large density noise and small positional noise means that the second term C In order to predict the scaling with system size of the smectic order parameters we determine in our simulations, we need to calculate |u(r+r , t)−u(r , t)| 2 ≡ |∆u(r)| 2 . This in general gets two contributions: a "bulk" term given by: and a "zero mode" term |∆u(r)| 2 0 given by: where the sum is, as usual, over all values q x = 2πm L , m = integer, allowed by our periodic boundary conditions, excluding m = 0.
The bulk term represents the contribution from all Fourier modes with q z = 0, while the "zero mode" term, as its name suggests, incorporates the contribution from modes with q z = 0. The latter modes would be absent in a system with fixed boundary conditions (i.e., u = 0 on the boundaries), but is present in our simulations, since we use periodic boundary conditions.
We'll consider first the contribution of the first term C ET (1) uu to |∆u(r)| 2 , which we'll call |∆u(r)| 2 1 , and then that of the second term C ET (2) uu , which we'll equally unimaginatively call |∆u(r)| 2 2 . Let's first consider the case L L c . In this case we are in the regime x L c and z L 2 c λ . In this regime, the range of wavevectors q x q c , q z q x dominates the contribution of the first term to both the bulk term (A.26) and the zero mode piece (A.27). We can therefore, for this range of r, replace the first term in (A.21) with Eq. (A.25), which, as noted earlier, has exactly the same form as in an equilibrium smectic, for which these correlations were computed long ago [25]. Hence, for the contribution from the bulk term from C ET (1) uu (q), we can simply quote the results of [25] with the trivial replacement of k B T in the equilibrium problem with ∆ u /2; this gives where we have defined the correlation lengths along x and z respectively, and the scaling function g(w) as a function of its argument w ≡ x 2 λ|z| is [25] The lengths ξ B x and ξ z are the distances along x and z at which the rms fluctuations |∆u(r)| 2 1,B are of order a l (more precisely, they are In most systems, the contribution of the zero modes are negligible, due to the factor of 1 L 2 in front of Eq. (A.27). In active smectics, however, if the system size L L c , the critical size above which active tension effects become important, so that we can drop the active tension term D ux q 2 x relative to the Kq 4 x term (since, for such system sizes, q x q c ≡ D ux /K), C uu (q x , q z = 0, t) is so large at small q x (diverging like 1 q 4 x as q x → 0), that this "zero mode" term can actually dominate. Indeed, for L L c , we find that, for x L, the sum in (A.27) is dominated by the smallest allowed q x 's, for which we can expand the cosine to leading order in x L . This gives where we have defined a correlation length ξ x given by such that when x = ξ x the zero mode real space fluctuations |∆u(r)| 2 0 = 2 , where q 0 ≡ 2π a l , with a l the layer spacing. This correlation length is clearly much shorter than that coming from the "bulk modes", which is given by Eq. (A.29), in the small noise (∆ u → 0) limit, since it diverges like 1/ √ ∆ u as ∆ u → 0, while the bulk correlation length diverges like 1/∆ u as ∆ u → 0. Hence, we expect ξ x to give the correlation length in the x direction for small displacement noise ∆ u . , which is clearly much smaller, for |z| 36λ, than the value of x at which the bulk term crosses over from being controlled by z to being controlled by x; i.e., |x| ∼ λ|z|. Hence, for all large r's at which the bulk term makes an appreciable contribution to |∆u(r)| 2 1 , we can replace that bulk term by 2 q 2 0 |z| ξz . Hence, we can always (for L L c ) and large r replace |∆u(r)| 2 1 = |∆u(r)| 2 1,B + |∆u(r)| 2 1,0 with Now we consider the case L L c . In this case the zero modes' contribution to |∆u(r)| 2 is always negligible. For x L c and z It is the integration of C ET (1) uu over small q's (i.e., q x q c ) that dominates in this regime of L. In this regime of q',s in the stable case D ux > 0, Γ u scales like q 2 for all directions of q. Specifically, in this limit, Γ u ≈ B q 2 z + D ux q 2 x . Using this in (A.23), and using the result in (A.26) plus the constant then gives where R ≡ x 2 + Dux B z 2 . Now let's consider the contribution of the second term in (A.21) to the real space fluctuations of u. Since this scales like 1/q 2 for all directions of q, it will also make contributions to the mean squared real space fluctuations |∆u(r)| 2 that scale like |∆u(r)| 2 2 ∝ ln r a l . This scaling holds for this term for all r. The detailed calculation goes as follows: we start by rewriting C ET (2) uu (q) in polar coordinates (q, θ) where (q x , q z ) = q(sin θ, cos θ): where we have defined: where Λ is an ultraviolet cutoff of order 1 a l . The integral over q in this expression can readily be evaluated for large r by dividing it, like Gaul § , into three parts: where we have defined II ≡ In these expressions, φ is the angle between r and the z-axis, and we have exploited the evenness of the cosine to replace its argument with its absolute value.
The integral II clearly converges as q → ∞, due to the oscillation of the cosine, combined with the 1 q falloff. We can therefore safely take the upper limit on this integral to ∞ for large r, where the lower limit is small, and obtain II ≈ The integral for Υ can be evaluated by the simple trigonometric substitution v ≡ tan θ, which gives which can be straightforwardly evaluated by complex contour techniques, giving: One might think that this logarithmic divergence would, for r L c , be dominated by the stronger power law divergence of the contributions |∆u(r)| 2 1 coming from the first term, at least for large r. But if the density noises ∆ ρ(x,z) are large and the positional noise ∆ u is small, this need not be the case, since the coefficient Υ of the logarithmic divergence has pieces that grow linearly with ∆ ρx,z . And we do in fact expect that, in our simulations, the density noises ∆ ρ are large while the displacement noise ∆ u is small, as explained in the main text.
So, in the limit of small angular noise, we expect ∆ ρ(x,z) to be large, by some suitable dimensionless measure, compared to ∆ u . This means that the contribution from the non-equilibrium second term in (A.20) to |∆u(r)| 2 can actually dominate that of the first term all the way out to the correlation length (i.e., ξ x given by (A.32) for L L c , or ξ B x given by (A.29) for L L c ) provided that the value of |∆u(r)| 2 2 is greater than 2 q 2 0 (the value of the contribution of the first term at the correlation length) there. This leads to the condition where X stands for either ξ x or ξ B x . Assuming, as seems reasonable, that the diffusion constantsD x,z and D ρ(x,z) are all comparable, and are larger than or comparable to D ux , and likewise assuming that B ∼ B and ∆ ρx ∼ ∆ ρz ≡ ∆ ρ , we can estimate Υ ∼ C∆ρ B D 2 , where D is the common value of the diffusion constantsD x,z and D ρ(x,z) . On dimensional grounds, we expect D ∼ ρ 0 C, where ρ 0 is the mean density. Using this estimate, we obtain a lower bound on the density noise required for the second, non-equilibrium term in (A.20) to dominate all the way out to the correlation lengths ξ x and ξ B x : which should be satisfied in our simulations since ∆ ρ is large, for the reasons discussed earlier, while ∆ u is small, which makes both ξ x and ξ B x large, thereby reducing the right hand side of (A.59).
In summary, for values of ∆ ρ that are large enough, in the sense just described, and values of the active tension D ux that are small enough, there are four important crossover lengths. The first is L c , beyond which the effects of the active tension become important. The other three, which are independent of D ux in this limit, are the correlation lengths ξ x , ξ B x (for L L c and L L c , respectively), and ξ z that give the distances along x and z beyond which |∆u(r)| 2 grows algebraically with x and z. These length scales divide the D ux -system-size L plane into five regions, as illustrated in Fig. 5. The active tension D ux is only important for length scales L L c ≡ K/|D ux |, which defines the upper and lower bounding curves. Because of this, regions I, II, and III extend right across the L-axis, where D ux = 0. In region I (L ξ x ), the second, entirely nonequilibrium, logarithmic contribution |∆u(r)| 2 2 to |∆u(r)| 2 dominates. In region II (ξ x L ξ z ), |∆u(r)| 2 is much bigger than the lattice spacing a l for r primarily alongx (i.e., the layer direction) and x ξ x , but is much less than a l for all points whose separation lies primarily in theẑ direction (i.e., the layer normal). Region IV (D ux > 0 and L L c ) is not further divided into more sub-regions by the correlation lengths ξ x and ξ z , since the behavior of |∆u(r)| 2 is always logarithmic for large r (i.e., x L c or z L 2 c /λ) in all directions, which solely determines the L-dependence of the quantitiy S n . And finally, for D ux < 0 and L L c , the smectic state is unstable. We do not have any analytic theory for the long-term, large distance state of the system in this regime. However, as argued in the main text, we expect the instability to induce nucleation of dislocations, which will cut off order prameter correlations at L c , hence, for L L c , we expect S n ∝ 1 A = 1 L 2 . To summarize the above rather complicated discussion, |∆u(r)| 2 = |∆u(r)| 2 1 + |∆u(r)| 2 2 , where |∆u(r)| 2 1,2 represent respectively the contributions of the first and second terms in (A.20). For L L c , |∆u(r)| 2 1 is given by while for L L c , it is given by And |∆u(r)| 2 2 is always given by regardless of whether L L c or L L c . We remind the reader that φ is the angle between r and the z-axis, Υ is a linear function of the noises ∆ ρ(x,z) and ∆ u given by Eq. (A.57), and the correlation lengths are given by ξ Appendix A.4. Smectic Order parameters, apolar case We'll now discuss the implications of these results for the smectic order parameters that we measure in our simulations. As in equilibrium, smectic order is characterized by an infinite set of local complex smectic order parameters ψ n (r, t) defined via [34] ρ(r, t) ≡ ρ 0 + ∞ n=1 ψ n (r , t)e inq 0 z + c.c. , where ρ is the number density, ρ 0 its mean, and q 0 ≡ 2π/a l , with a l the distance between neighboring layers. To quantify the order in our simulations, we numerically determine the quantities S n ≡ I(q n , t), where the "intensity" ‡ I(q, t) ≡ |ρ(q,t)| 2 N 2 , with N = ρ 0 A the number of particles in the system, A = L 2 is the area of the system, q n ≡ nq 0ẑ , and ρ(q, t) = d 2 r ρ(r, t)e −iq·r . (A.64) In writing the expansion (A.63), we implicitly assume that the ψ n (r, t) are slowly varying in space and time (i.e., that they only have support in Fourier space at small wavevectors and frequency). That is, we have absorbed the rapid variation of the density in the smectic phase into the complex exponentials e inq 0 z , while the ψ n (r, t) embody the ‡ So-called, because it is proportional to the intensity of scattering of electromagnetic radiation by the smectic when the wavelength of that radiation is comparable to the smectic layer spacing. slow spatial variations of the positions u(r, t) of the layers. Given this, the definition (A.63) is equivalent to: from which it follows from our definition of the S n 's that The phases of these local order parameters are all proportional to the displacement field u; that is, we can write ψ n (r, t) = |ψ n (r, t)|e −inq 0 u(r ,t) . Hence, since the amplitudes |ψ n (r, t)| of the smectic order parameters are not Goldstone modes of the system, and are therefore not expected to have large fluctuations, the decay of the of the correlation function ψ * n (r, t)ψ n (0, t) with r is driven primarily by the fluctuations of the layer displacement u(r, t). Motivated by this, we replace |ψ n (r, t)| with a real positive constant ψ 0 n , and obtain for the correlation function we need: This can be related to the mean squared real space fluctuations |∆u(r)| 2 by noting that ∆u(r) is a zero-mean Gaussian random variable. This follows from the fact that it is a linear function of the Fourier components u(q, ω), which are in turn linear functions of the Fourier components f u (q, ω) and f ρ (q, t). Since these noises are themselves, by assumption, zero-mean Gaussian random variables, so is ∆u(r − r ).
Using this fact, we can use the well-known (and easily derived by any who don't know it) relation for any zero-mean Gaussian random variable x that e ix = exp(− For L ξ x and L L c , denoted as region I in Fig. 5, the exponential factor in this expression is nearly 1 over the entire region of integration (recall that ξ z ξ x ), and, as a result, the S n 's fall off algebraically with L: S n ∝ L −n 2 η I . For ξ z L ξ x and L L c , denoted as region II in Fig. 5, we can ignore the z-dependence of the exponential factor, but not the x-dependence. The integral is then clearly dominated by z ∼ L ξ x and x ∼ ξ x , which means x z in this dominant region. We can therefore replace r with z and φ with 0 in (A.71), and obtain Finally, for L ξ z , the integral over x and z in this expression converges as L → ∞; hence, in this regime, denoted region III in Fig. 5, all of the S n fall off with system size like 1 A = L −2 . (Note that this regime must also have L L c ). For L L c in the stable regime D ux > 0, denoted as region IV in Fig. 5, both contributions |∆u(r)| 2 1 and |∆u(r)| 2 2 to |∆u(r)| 2 grow like ln r a l for large r, and hence, like region I, this region will also exhibit algebraic decay of the S n 's: S n ∝ L −n 2 η IV , albeit with η IV now taking on a different value: However, since ∆ u is expected to be small, as discussed in Sec. Appendix A.2, the numerical values of η IV and η I are very close to each other, as shown in the simulation. For L L c in the unstable regime D ux < 0, denoted as region V in Fig. 5, the system is unstable, and we expect dislocations to proliferate, cutting off order parameter correlations at L c , and causing S n to again decay like L −2 .
Appendix B. Hydrodynamic theory predictions for the polar active smectic P phase Appendix B.1. Linearized Eigenfrequencies and Instability threshold Fourier transforming Eqs. (3.9),(3.10) in the main text and only keeping terms to linear order in u and δρ, we obtain with Γ u (q) and Γ ρ (q) having exactly the same expressions as in the apolar case. Guessing that one of the eigenmodes has eigenfrequency with α(q) = O(q 2 ), and inserting this guess into Eq. (B.1) with f u set to zero on the right hand side, we see that, to leading order in q, u = Cqz vρqx δρ. Inserting this into equation (B.2) with f ρ set to zero on the right hand side gives where we have defined Thus, Now looking for a second mode with eigenfrequency ω 2 = O(q 2 ), Eq. (B. 2) with f ρ set to zero on the right hand side, we see that, to leading order in q, δρ = − iv 2 vρ q z u. Inserting this into equation (B.1) gives where we have defined For q x q c , where the D ux term is negligible relative to the K term, Γ u2 (q) again looks like an inverse smectic propagator, scaling like 1/q 2 z for q z λq 2 x , and like 1/q 4 x for q z λq 2 x . Note that only ω 2 involves the active tension D ux . Hence, only this eigenfrequency can acquire a positive imaginary part, signaling an instability, when the active tension D ux goes negative. And since the imaginary part of ω 2 is manifestly an increasing function of q 2 z , it is obvious that this instability must first set in at q z = 0, q x = q c = −D ux /K.

Appendix B.2. Linearized u − u correlation functions
With these eigenfrequencies in hand, we can now compute the spatiotemporally Fourier transformed u − u correlation function in the linearized approximation precisely as we did in the apolar case. This gives This can again be integrated by parts to obtain the equal-time u−u correlation function, albeit not as neatly as in the apolar case. We find, after considerably more algebra, Note that once again, as in the apolar case, the first term will make a contribution to |u(r)| 2 that is ∝ √ L (for L L c ). Unlike the apolar case, however, the second term now only makes a finite contribution to |u(r)| 2 . This can be seen by noting that for most directions of wavevector q, the second term is independent of the magnitude q of q, since both the numerator and denominator of the second term scale like q 6 for such generic directions of q. This scaling does not hold, however, for q x q 2 z . In this regime, the second term is readily seen to be ∝ 1/q 2 . However, there is not enough phase space in the regime q x q 2 z for this divergence to lead to any divergent fluctuations in real space; that is, does not diverge as L → ∞.
The convergence of this second term means that its contribution to the real space u − u correlations can be dropped. We are thus left with only the first term in (B.10), which has exactly the same structure as the corresponding term in the apolar case. A moment's reflection reveals that this must mean that all of the scaling regimes found in the apolar case also exist in the polar case, in the linear approximation to the latter. The only difference is in region I, where the algebraic decay of order that was induced by the second term in the apolar case is absent in the polar one.

Appendix C. Nonlinear effects for polar smectic P phases
Going beyond the linearized approximation, the polar case has symmetry-allowed nonlinear terms in its equation of motion (i.e., those explicitly displayed in the equations of motion (3.9), (3.10) in the main text). These become important only at the extremely large length scale L N L discussed in the main text, and lead to the appearance of region VI in Fig. 5 for the polar case. We have not been able to determine the behavior of the smectic order parameters in this regime, but strongly suspect, as discussed in the main text, that they fall off as L −2 with increasing system size L in this regime, which includes the L → ∞ limit for the stable case.
The discussion in Secs. Appendix B.1 and Appendix B.2 was based entirely on the linear approximation to the full equations of motion (i.e., Eqs. (3.9), (3.10) in the main text). As mentioned earlier, in the polar case the non-linearities in these equations are relevant (technically, "marginal") in d = 2. These terms give rise to fluctuation-induced departures from the linearized theory that are proportional to the square of the noise strength σ (note that our hydrodynamic parameters ∆ u,ρx,ρz are all proportional to the square of the noise strength σ, as that noise strength is defined in equation (2.2). To see this, note that, since the speed of the particles is fixed, and the angular noise in (2.2) is proportional to σ, then the velocity noise f will be of order σv 0 , and, hence, also proportional to σ. Hence, the mean-squared forces, which are proportional to the various ∆'s, are proportional to σ 2 .) Furthermore, in d = 2, they grow logarithmically with system size L. Hence, for small noise strengths σ, the system must be larger than a nonlinear length scale L NL which grows very rapidly with decreasing noise strength: L NL ∼ exp (constant/σ 2 ) before the linear theory just described becomes invalid.
What happens beyond this long length scale L NL remains an open question, which can only be answered by a full renormalization group analysis. Due to the large number of marginal non-linearities (six) in our problem, such an analysis is quite formidable, and we have not attempted it. Here we will restrict ourselves to plausible speculation.
One possibility is that the non-linearities eventually induce the unbinding of dislocations, thereby making the smectic order parameter fall off far more rapidly than algebraically (presumably exponentially) for longer length scales. This would be consistent with our observation that, in the largest systems we study, dislocations do appear, and the smectic order parameter does fall off exponentially with increasing system size.
Independent of the question of dislocations, however, one thing that will certainly happen for L L N L is that the scaling law predicted by the linear theory relating higher harmonic smectic order parameters S n to the lowest harmonic order parameter S 1 , namely S n ∝ S n 2 1 , will break down, since that scaling law was derived in the linear approximation, in which the fluctuations of the layer displacement u(r, t) are Gaussian. Once non-linearities become important (which, by definition, happens for L L N L ), all of the fluctuations will become non-Gaussian, and this scaling law will break down.
Therefore, the hydrodynamic theory predicts that the S n ∝ S n 2 1 for L L N L , and that this law will break down for L L N L . Since, in light of above discussion, we also expect correlations to fall off more rapidly for L L N L , this means we expect both a breakdown of scaling law S n ∝ S n 2 1 relating the smectic order parameters, and the power law scaling S 1 (L) ∝ L −η I to fail at same system size L (namely, L ∼ L N L ). And this is exactly what we see in the simulations, as illustrated in Figure 1g of the main text, where we again plot S 1 (L), S 2 (L) 1/4 , and S 3 (L) 1/9 versus system size L. These lie on top of each other at small distances, consistent with the linear theory, but depart from each other at longer length scales, signaling the onset of important non-linear effects at large distances. And at the system size L where this happens, the order parameters (particularly S 1 ), begin to depart from power law scaling with system size L (see again Figure 1g of the main text). This is entirely consistent with our argument above about the impact of the relevant nonlinearities on the order parameter correlations. We therefore interpret the system size at which both these phenomena happen as L N L for our system.
It would be a very dramatic illustration of this fundamental difference between the apolar and the polar case to show that Eq. (3.8) in the main text holds out to arbitrarily long length scales in the apolar case, which the theory says it should, since there are no relevant nonlinearities allowed by symmetry in that case. However, this can only be tested in systems in which the apolar case is actually stable (so one can explore region IV in Figure 3). Alas, as discussed earlier, all of the apolar systems we have investigated appear to be unstable. For numerical convenience, we have accomplished this by adjusting the interaction range r int while keeping the system size fixed, until the average layer spacing a converged to unity. Only in the case of isotropic repulsion was this necessary, as only in that case was the average layer spacing significantly lower than the interaction range (specifically, we found a l ≈ 0.95 r int ). The modified values for the interaction range needed to increase a l back to 1 for isotropic repulsion were r int = 1.05 for the F-model and r int = 1.06 for the N-model. In order to keep the speed of individual particles comparable across different repulsion types, the speed (displacement per time step) was always measured in units of the interaction range. In simulations, this choice of a l = 1 made it desirable, for the reasons just discussed, to always work with systems of integer size L. This not only minimized potential mismatch between the layer spacing and the system size, as discussed above, but also ensured that the spatial subdivision algorithm employed by us for efficient calculation of local interactions in large systems functioned correctly.

Appendix G. Winding number
Periodic boundary conditions make layered patterns with different integer "winding numbers" w possible: w = 0 means each layer connects to itself at the periodic boundary, while for w = 1 each layer connects to the one above it, w = 2 implies connecting to the layer two above it, etc. To minimize this artifact of periodic boundary conditions, we restrict our simulations of smectic patterns to "flat" (i.e., w = 0) configurations. However, it should be noted that in the apolar system, the undulation instability at large length scales breaks the smectic layers, which may lead to the emergence of a smectic pattern with different winding numbers.