A growing bacterial colony in two dimensions as an active nematic

How a single bacterium becomes a colony of many thousand cells is important in biomedicine and food safety. Much is known about the molecular and genetic bases of this process, but less about the underlying physical mechanisms. Here we study the growth of single-layer micro-colonies of rod-shaped Escherichia coli bacteria confined to just under the surface of soft agarose by a glass slide. Analysing this system as a liquid crystal, we find that growth-induced activity fragments the colony into microdomains of well-defined size, whilst the associated flow orients it tangentially at the boundary. Topological defect pairs with charges \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm {\textstyle{1 \over 2}}$$\end{document}±12 are produced at a constant rate, with the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+ {\textstyle{1 \over 2}}$$\end{document}+12 defects being propelled to the periphery. Theoretical modelling suggests that these phenomena have different physical origins from similar observations in other extensile active nematics, and a growing bacterial colony belongs to a new universality class, with features reminiscent of the expanding universe.

A ctive matter occurs at many length-scales: from bird flocks 1 , through shaken grains 2 and swimming bacteria 3 and synthetic colloids 4 , to gels in which protein motors 'walk' on filamentous 'rails' 5 . All contain agents consuming energy to drive locomotion. The field attracts physicists (as a 'grand challenge' in non-equilibrium statistical mechanics), theoretical biologists (as a paradigm uniting phenomena across disparate scales) and materials scientists (who see applications in, e.g., self-assembly and artificial wound healing).
Active matter systems may show 'polar' (2π-periodic, like an arrow) or 'nematic' (π-periodic, like a double arrow) orientational ordering. Within the part of the system under consideration, momentum is either conserved or not. A bulk system which conserves momentum is termed 'wet', while in a 'dry' system, such as a two-dimensional (2D) layer interacting strongly with a substrate, momentum is not conserved. The polar/nematic and wet/dry dichotomies make up the four classical universality classes of active matter 6 . Recent developments show that 'wet' and 'dry' are the limiting cases of a continuum of possible theories 7 , and that 'scalar active matter' exists with no orientational order 8 . In all these cases, the number of active agents is conserved at all times. Studies of biological morphogenesis as active matter have also, until recently 9,10 , focussed on developmental phases with constant average particle number by balancing cell birth and death. In the following, we refer to these systems generically as 'numberconserving' active matter.
Here we instead study a simple active matter system with nonconserved particle number: a two-dimensional layer of dividing Escherichia coli cells confined in agarose gel by a glass slide. It represents the early stage of bacterial colonisation in many contexts, and constitutes perhaps the simplest example of biological morphogenesis. Despite the absence of cell differentiation, a growing bacterial colony is governed by similar physical constraints to those controlling eukaryotic morphogenesis, including tumour growth 11 .
Cell-cell mechanical interactions constrain colony growth [12][13][14][15][16][17][18][19] . An elongating E. coli pushing against its neighbours in a colony has been modelled as an extensile dipole in a viscous liquid, i.e., a wet active nematic, in which viscous drag between the cells and bulk liquid controls colony morphology 10 . A dry continuum model has also been proposed to describe the chaotic behaviour of nematic domains in simulated colonies 19 . However, no theory of growing active nematics has yet been tested against quantitative experimental data.
We provide such data by analysing growing E. coli colonies as 'living anisotropic fluids'. We observe phenomena resembling those seen in number-conserving active nematics, including 'active anchoring' at the boundary 20 and the proliferation of topological defects 21 . However, the underlying mechanisms differ. We propose a theoretical framework in which our observations emerge from growth-driven Hubblelike expansive flows 22 and friction between cells and their substrate 17 . Alluding to this qualitative broad analogy, henceforth we call this class of systems (dry) 'Hubble active nematics'. [Clearly, this suggestive analogy comes with limitations. Most notably, within our quasi-2D colonies the amount of matter grows over time, where in the expanding universe it is the metric that changes in scale.]

Results
Morphology and growth rate. We begin by reporting results on the morphology of our growing microcolonies (see Methods for experimental details on growth conditions and image processing). Confocal imaging ( Supplementary Fig. 1) shows that cells are embedded just underneath the agarose surface, confined on one side by the cover slip. Figure 1 shows snapshots from a typical colony, which are representative of the behaviour we observe under our experimental conditions (see Methods). A number of qualitative changes occur as cells lengthen and divide. At early times, the colony is elongated, Fig. 1a, b, as cells grow longitudinally along a common axis. Later, this 1D order breaks down, giving way to a pattern of nematic domains bordered by defects, Fig. 1c. During this stage, cells push outwards in all directions, and the colony becomes more isotropic, Fig. 1 (inset, see also Supplementary Note 1 and Supplementary Fig. 2). Under our conditions, the middle of the colony turns phase dark after~5 h when it contains N ≲ 10 3 cells, Fig. 1d (dark spots), indicating buckling of the 2D layer into the third dimension 17 , and we stop our analysis. Typically, buckling occurs when the colony size 2R ≳ 50 μm. All colonies show the same trend, although quantitative details are subject to noise ( Supplementary Fig. 2).
Cells in the middle of a large colony may grow slower due to nutrient limitation. Measurement of cell lengths as a function of time, Fig. 2a, shows that their growth rate is in fact spatially uniform to within experimental uncertainties, so that nutrient limitation is unimportant for our colonies (although the growth rate is time dependent, Supplementary Note 2 and Supplementary Fig. 2b). The near-uniform growth rate is also the reason why we do not observe fingering instabilities, which appear when growth rate is limited by nutrient 16 or depends on local density 23 .
'Hubble' expansion. We found no coherent azimuthal movement of cells, but there is a coherent radial movement, whose speed increases linearly with distance from the centroid, Fig. 2b, reminiscent of the Hubble expansion of the universe. Consider the velocity field, v, of this growth-induced motion: we now write an equation for this quantity based on our experimental observations. Mass conservation requires that the cell density, ρ, satisfies where Λ is the growth rate, and D t = ∂ t + v · ∇ is the material derivative. Assuming incompressibility (we relax this assumption below), i.e., D t ρ = 0, we find where hereafter repeated suffices are summed. Since Λ is independent of space in our colonies, and there is no coherent azimuthal movement, Eq. (2) predicts the radial component to be v r = Λr/2 ≡ Hr, as observed, Fig. 2b (the factor 1/2 arises as the system is 2D). A linear fit to the population-averaged data in Fig. 2b gives our 'Hubble constant' H = 0.007 min −1 . Theoretically, H = 0.5Λ, while individual values of (H, Λ) obtained from each of the 32 colonies, Fig. 2c, are consistent with a linear scaling of H with Λ and are fitted by H = 0.4Λ (the correlation between H and Λ is r 0.77, corresponding to a p-value p < 0.00001). The discrepancy from the theoretical prediction H = 0.5Λ may arise either from small systematic errors in the fitting (e.g., due to small drift close to colony centroids) or from the approximations in the theoretical argument, such as that of isotropic growth. The scatter in the data in Fig. 2c gives an indication of the extent of intrinsic noise in our population of colonies.
Note that the expansive 'Hubble flow' we observe is distinct from the spontaneous shear flow exhibited by number-conserving active nematics 24,25 .
Active anchoring. A striking feature of the larger colonies is that cells tend to be aligned tangentially at the periphery, Fig. 1d. We fitted a Bézier curve through the centroids of the outermost cells, and measured the angle ψ between the curve and each boundary cell. For a typical colony just before it buckled, most peripheral cells show ψ ≲ 30°, Fig. 3a. Such 'active anchoring' is also seen in simulated number-conserving active nematics 20 , where it originates from the ingress of topological defects following the buckling of the boundary between an active nematic and an isotropic continuum. We will see later that the mechanism in our case is different.
Continuum description. Turning to the bulk of our colonies, we first construct a continuum structural description based on our raw data, which consist of the position, length and orientation of each cell as functions of time, {r i (t), L i (t), ν i (t)}. We start this coarse graining procedure by defining a function that 'smears out' individual rod-shaped bacteria of diameter W: with Δr i = r − r i , v i the unit vector perpendicular to ν i and σ a smoothing length. We use σ = W to probe ordering at the singlecell level (σ → 0 gives sharp L i × W rectangles). The integrated cell 'mass' is preserved as σ varies. [Note that the functional form of f i we use is required to account for the non-spherical, rod-like shape of bacteria.] The (number) density field is then  To quantify orientational order, we assume that the two poles of a cell are indistinguishable (i.e., ±ν i are equivalent), and construct a traceless and symmetric Q-tensor 26 field where we have also defined the scalar order parameter and the director n(r), which is the unit eigenvector associated with the positive eigenvalue of Q αβ (r). Supplementary Fig. 3 (see Supplementary Note 4) shows Q αβ (r) of the snapshots in Fig. 1.
Global order. The global version of Eq. (7) quantifies the degree of orientational order averaged over an N-cell colony: where θ i is the angle between the axis of cell i and the director (i.e., ν i × n = cosθ i ), and 〈…〉 denotes averaging over a colony. This evaluates to unity when all cells are perfectly aligned, and to 1= ffiffiffiffi N p for total orientational disorder. The measured degree of global order, S(N), varies significantly between colonies, Fig. 3b, but the trend is clear. Global order is high for up to ≳3 generations (N ≲ 10), and then decreases. Significantly, forms the lower envelope of the vast majority of our data. Thus, while there is no sustained growth-induced global ordering 14 , substantial residual order persists up to the buckling point, when our observations stop.
Correlation lengths and domains. Underlying the persistence of residual global order is a significant degree of local orientational order in the form of domains of aligned cells: Fig. 1c, d. Such local order can be quantified by calculating two-point correlation functions of the Q-tensor, C || (r) and C ⊥ (r) (Supplementary Note 5), which measure the azimuthally-averaged orientational correlation at a distance r from a typical cell parallel and perpendicular to its axis, respectively. These functions for the colony snapshots in Fig. 1 are shown in Supplementary Fig. 4. In the earliest stage of growth (N = 15), C || (r) and C ⊥ (r) decay over ≈10 and 5 μm respectively, i.e., the whole colony is ordered (as is obvious from Fig. 1a). At intermediate times, parallel correlations grow faster than perpendicular correlations, before C || (r) and C ⊥ (r) become more comparable at late times. This trend is generic. We calculated C || (r) and C ⊥ (r) for 32 colonies, and obtained correlation lengths defined by l jj;?
so that two cells each of length and width L ≳ 2 μm and W ≲ 1 μm placed side by side and end-on give l jj c ¼ L=2 % 1 μm and l ? c ¼ W=2 % 0:5 μm respectively. At early times, l jj c =l ? c grows ( Supplementary Fig. 5), with l jj c growing faster and peaking at ≲5 μm, compared to a maximum l ? c ≲ 3 μm, corresponding to domains of some tens of cells, Fig. 1b. Thereafter, both correlation lengths decrease to ≳2 μm, corresponding to domains of ≲10 cells. Similarly-sized 'micro-domains' have been seen in recent simulated 2D bacterial colonies at late times 19 .
To highlight the role of active growth in generating this pattern of behaviour, we performed Monte Carlo simulations [27][28][29][30][31] of passive (non-growing) sphero-cylinders 32 with lengths and widths taken from bacterial colonies (see Supplementary Note 6 and Supplementary Fig. 6). In the passive case, l jj c =l ? c % 1 throughout the entire history of the colony ( Supplementary  Fig. 5).
Proliferation of topological defects. The existence of domains implies topological defects, where the director field is illdefined. Defects control the elasticity and dynamics of passive liquid crystals, and feature prominently in numberconserving active nematics 6 , where they arise from the self-sustaining shear flow arising from the stirring of the fluid by active particles.
The topological charge, or winding number, of a defect, is the accumulated change (in units of 2π) in the angle ϕ made by the director field relative to an arbitrary direction on following a closed loop, C, around the defect: Defects are identified via a topological charge density 20 : which can be computed directly from the Q-tensor field. An example of spatial maps of q is shown in Supplementary  Fig. 7 (see also Supplementary Note 7). We identified individual defects as local maxima and minima of q(r), and determined their polarity from the unit vectorp d;α ¼ À∂ β Q αβ = ∂ β Q αβ at the local maximum. As in conventional active nematics, defects in our colonies mostly have topological charge ± 1 2 . As the colony grows, we see the formation of multiple defect pairs in the colony interior, Fig. 4a. While defects are created continuously, the total signed topological charge R qd 2 r ¼ 0. However, the total number of defects, which we compute as the total unsigned topological charge, R q j jd 2 r scales near-linearly with N, Fig. 4b, suggesting that defect generation and annihilation is a bulk, not boundary, phenomenon, and that topological (unsigned) charge grows exponentially in time. Linear scaling with N is expected on general grounds if, as we observe, the colony fragments into micro-domain of well-defined size (≲10 cells).
Defect dynamics. Tracking the topological defects through time, Fig. 4a (see also Supplementary Fig. 7), shows that the orientation of the three-fold symmetric À 1 2 defects does not significantly impact their dynamics 21 . Presumably, these are simply advected by the expansion of the colony, perhaps aided by mutual elastic interactions. In contrast, comet-like 'head and tail' þ 1 2 defects show marked directed motion-the speed just prior to buckling is v +~0 .05 μm s −1 . [This value should be seen as an order-ofmagnitude estimate as the measurement requires unambiguous tracking of defects which we could only do in few selected cases.] The distribution of b p d Á Δb r dÀcom , Fig. 4c, where Δb r dÀcom is the unit vector pointing from a colony's centre of mass to the defect, shows that positive defects propel along the direction of b p d (i.e., tail to head), reminiscent of their behaviour in a non-growing extensile active nematic 21 .
The spatial distribution of ± 1 2 defects, Fig. 4d, shows 'charge separation'. There is an excess of À 1 2 defects near the colony centre, while þ 1 2 defects predominate at the periphery. Presumably, inwardly-directed þ 1 2 defects have a high chance of annihilating À 1 2 defects advected by colony expansion, while outward-directed þ 1 2 defects rapidly propel to the colony edge, where advected À 1 2 defects have not yet reached. Moreover, once a þ 1 2 defect has propelled itself to the colony edge in a tail-to-head direction (Fig. 4c), the director field pattern near its 'head' matches the tangential anchoring, Fig. 3a, and is therefore energetically favoured. Thus, we find that þ 1 2 defects are rarely ejected from the boundary into the colony. This contrasts with a number-conserving active nematic surrounded by an isotropic medium, where ejection of þ 1 2 defects from a deeply-undulating boundary into the interior is prominent 20 .
Intriguingly, the region in which we observe an excess of À 1 2 defects, r ≲ 0.2R (here R is the radius of the inscribed circle), is precisely where we see the onset of buckling into the third dimension. This is reminiscent of the physics of passive liquid crystals, where the director field can avoid the formation of a costly defect core by escaping into the third dimension. This observation also suggests that defects may provide a means to relax stress in the system, as proposed recently for eukaryotic colonies 33 .
Theoretical description and length scales. Our observations and analyses so far suggest that it should be fruitful to describe a growing 2D colony of E. coli cells theoretically as an active nematic, albeit in a different universality class than almost all current theories of this kind because of particle number nonconservation. Here, we first delineate a number of components that must feature in such a theory and point out their implications. We then perform numerical simulations of a set of equations which assumes a particularly simple constitutive relation for the active stress. There are two important length scales in our system. The first one, l a , is related to the onset of the 'generic instability' in active nematics.In wet systems, l a is given by a balance between active and elastic stresses 34 , as l a $ ffiffiffiffiffiffiffiffi ffi K=a p , where K is a Frank elastic constant and a is the active stress coefficient, defined so that the growth-dependent active stress scales as aQ αβ . In a growing colony, K and a both depend on growthmediated cell deformations, and should scale as the cell wall elastic modulus, G, so that a~G and K~GL 2 for cells of length L. They should also scale similarly with cell density 16,19 , so that we may expect K/a to be a constant, and l a~L . In our dry system, the functional dependence of l a on parameters is more complicated (see below and 35 ). However, as l a is generally lower for the dry than for the wet limit 7 , we expect l a~L should still hold for our colonies.
A growing cell pushes and deforms its neighbours, generating stresses that scale as G. The resulting cell movements are opposed by dissipative forces, which we take to be lubricated frictional sliding between cells and agarose or glass controlled by a frictional drag coefficient per unit length, γ, so that the drag force on a cell of length L sliding at speed v is~γLv. Force balance in terms of the stress tensor σ αβ then reads where ρ is the cell density. Our 'dry' description in Eq. (11) is motivated by the fact that our colony is very close to the glass coverslip ( Supplementary Fig. 1); in general this equation will hold above a hydrodynamic 'screening length' l s $ ffiffi η γ q , with η the effective viscosity of the colony.
Growth stress propagates by diffusion; dimensional analysis gives the diffusivity $ G=γ 16 . [More rigorously, we can show that, within our theory, when Λ = 0 both ρ and the stress evolve according to a diffusion equation where the diffusion coefficient is $ G=γ, given the simple constitutive relation for the stress considered below.] The extent of stress propagation during the lifetime of a single cell (birth to division) defines a stress propagation length Previous measurements give G~10 6 Pa, albeit with large variations, while we estimateγ $ 10 15 (Pa s)m −2 using previous experimental results 36

(Supplementary Note 8). Thus l σm
m. This is an upper bound for the size of a colony in which internal growth stresses can relax into the agarose matrix. In practice, buckling into the third dimension will occur when the stress approaches the yield stress of 2% agarose, $ 3 10 3 Pa ( G. With this revised stress estimate, l σ~5 0 μm, similar to experiments. Intriguingly, v + Λ −1~5 0 μm as well, suggesting that defects may have a role in relaxing growth stresses. Growth-induced alignment. To see how growth should induce orientational order, consider first a system of bacteria oriented along y (Q yy = 1, Q xx = −1, Q xy = Q yx = 0), Fig. 5a, with a single slightly-misaligned bacterium B (centroid at y 0 ) whose angle θ with the orthogonal (x) axis is <90°. The only velocity gradient component is ∂ y v y = Λ. The fraction of B at larger y moves with higher v y , and so must experience a higher force than the lower part; the net torque realigns this bacterium.
Quantitatively, the force giving rise to the velocity of a segment, dl, of B located a distance l from the centroid must balance the frictional drag with the substrate at this position, γv y (y)dl. The total torque on B is then Setting τ = ζdθ/dt, where ζ is the rotational drag on a single cell due to its neighbours, and expanding to first order gives where ∂ y v y = Λ > 0, Eq. (14) has a stable equilibrium at θ = 90°, i.e., growth drives alignment. This quasi-1D analysis correctly predicts a high degree of alignment in our early-stage colonies, Fig. 1a, b. This persists only for a few generations, because the active instability length scale l ã L. When the colony size 2R increases to the point when 2R=l a ) 1, the extensile stresses due to growth create bend distortions, as in number-conserving active nematics. The colony now grows in two dimensions.
For simplicity, consider the simplest case of circularly symmetric expansion in 2D, for which the cell velocity field is v ¼ gðrÞr; with g(r) a local growth rate. Substituting Eq. (15) into nematohydrodynamic equations describing flow-orientation coupling (see Supplementary Note 9) yields the analogue of Eq. (14) for 2D: where ϕ denotes the azimuthal angle in polar coordinates and ξ > 0 is a flow-aligning parameter. Equation (16) predicts that growth-induced alignment depends on the sign of g′. For compressible flows, g′ > 0, there is a stable steady state (∂ t θ = 0) with θ = ϕ, giving a radially-ordered colony with a central +1 vortex defect. This is the analogue of the growth-induced alignment predicted by Eq. (14) for 1D.
We do not observe such alignment, because experimentally, we find incompressibility, i.e., g ≃ Λ/2 and g′ = 0, Fig. 2a. In this case, the steady state of Eq. (16) does not require any particular relation between θ and ϕ. The colony breaks into domains with a length scale controlled by l a . The number of defects scales as the number of domains. There are~10 ≈ 3 2 cells per defect, Fig. 4b, consistent with l a~L . Our colonies are effectively incompressible because their sizes (2R ≲ 50 μm) are always < l σ , so that growthinduced stresses relax quickly.
At the edge of a colony, g ≃ Λ/2 (inside) transitions to g = 0 (outside), so that g′ < 0. Equation (16) predicts a stable steady state with θ = ϕ + π/2, i.e., tangential alignment, again as observed. Qualitatively, the origins of this effect is easy to understand. The presence and absence of growth-generated forces inside and outside the colony respectively, Fig. 5b, means that a cell not aligned tangentially will be rotated back into such alignment, as previously suggested 37 .
Non-uniform active stress. In theories of number-conserving active nematics, the active stress is taken to be where a is a positive or negative constant for extensile or contractile activity respectively 6 . It is tempting to model growing bacterial colonies in the same way with a ∝ Λ. To see why this does not work for a Hubble active nematic, consider again a 1D colony in the y direction. Force balance, Eq. For a 1D colony occupying 0 < y < y 0 with spatially-constant Λ (cf. Fig. 2a), Eq. (18) solves to 14,15 σ yy ¼ ÀP 0 À 1 2 Λγ y 2 0 À y 2 where P 0 is the pressure exerted by the agarose at the boundaries. The active stress field is non-uniform, increasing from the periphery to the centre by Δσ $γΛR 2 for a colony of dimension 2R. Numerically, this inhomogeneity can be substantial (our estimated parameters give Δσ~3 × 10 3 Pa for 2R~l σ ). Physically, such non-uniformity necessarily arises because peripheral cells have to push at fewer cells sliding frictionally against the substrate as the colony expands 17 . The coupling of mechanics and nematic order in 2D gives rise to new complexity. Inspired by Eq. (19), we may expect A simple constitutive relation consistent with a dry compressible active nematic is p = G max[(ρ/ρ 0 − 1), 0], and a = a 0 ρ, where ρ 0 is the close packed density of undeformed cells. The equations resulting from substituting this or other constitutive relations into the laws of nematohydrodynamics are only tractable by numerical simulations (see below). Nevertheless, the same physics we appealed to in the 1D case will still give rise to a stress field increasing from periphery to centre.
In a sufficiently small colony, 2R ( l σ , these non-uniform growth-induced stresses relax rapidly through the periphery into the agarose, so that the flow is effectively incompressible, Fig. 2b. As 2R → l σ , however, stress relaxation becomes progressively more incomplete, eventually leading to the buckling when the colony is not confined to 2 dimensions (as is the case in our experiments).
Our simulations confirm all qualitative predictions of our theory, and also compares favourably with experiments, see Fig. 6 and Supplementary Fig. 8. First, as expected the radial velocity increases linearly with distance to the colony centre-interestingly, our simulations also typically show a slighly smaller 'Hubble constant' than expected, as in experiments (H < Λ/2, Figs. 2b and 6a). Second, simulations reproduce the tangential alignment at the colony boundaries ( Fig. 6b and Supplementary Fig. 8). Third, when the active coefficient a 0 is sufficiently large, the colony breaks up into domains of well-defined lengthscale, l a , separated by defects ( Fig. 6b and Supplementary Fig. 8). Fourth, the total topological charge increases near-linearly with colony area (proportional to number of bacteria, see Fig. 6c). Finally, the positive topological charges (mainly associated with +1/2 defects) are more peripheral than the negative ones (mainly associated with −1/2 defects, Fig. 6d).
An analysis of our data suggest that there is a substantial range of a 0 for which l a $ a À1=2 0 ( Supplementary Fig. 9). This is in line with simulations of turbulent-like dynamics in numberconserving dry active nematics 35 . Comparison with the latter study and with the results of 19  . A more systematic study would be required to confirm this however. As a 0 decreases, the domain size increases until for sufficiently small a 0 the bending instability is no longer active and the colony resembles a nematic tactoidal droplet with tangential anchoring (Supplementary Fig. 8).

Discussion
We have described the structural evolution of a growing 2D colony of rod-shaped E. coli until the point it buckles into the third dimension using the language of liquid crystal physics. Our system displays features that appear similar to those that are familiar from number-conserving active extensile nematics, such as tangential anchoring at the periphery, Fig. 3a, and a proliferation of ± 1 2 defects, Fig. 4. However, the underlying physical mechanisms differ. Thus, e.g., while invaginations from the boundary is a dominant mechanism of defect generation in a number-conserving active nematic embedded in an isotropic medium, we find that defects proliferated from the bulk of our colonies. The expansive growth of cells in these colonies drives new physics; we have named this system a 'Hubble' active nematic.
We have also proposed a continuum theory of such systems. Within this framework, the physics is controlled by two length scales: l a associated with the generic bending instability in dry active nematics, and l σ associated with the diffusive propagation of growth stress. For a colony with size 2R ( l σ , growth stress relaxes rapidly, and the cell flow is effectively incompressible. The colony breaks up into domains of size~l a . A colony with 2R → l σ will accumulate stress and eventually buckle at the centre, where the accumulated stress is highest. In all cases, we expect tangential alignment at the periphery. All these predictions are in agreement with our experimental observations, and are confirmed by computer simulations of this model. More quantitatively, such simulations further reproduce our experimental finding that the total number of defects increases approximately linearly with colony area, and that positive charge defects are more peripheral than negative charge ones.
To enable even more quantitative comparison between experiments and simulations of our growing active nematic system, progress needs to be made on a number of experimental fronts. Perhaps most importantly, the nature of the friction between growing cells and their substrate, i.e. the physics underlying our parameterγ, needs to be elucidated. There is also a need to understand the mechanics of bacterial growth into bulk agarose. It is known that the elastic properties of the surroundings have a complex effect on, e.g., the point at which a 2D colony buckles into the third dimension 17 . The detailed mechanisms for these effects are far from understood. The effect of bacterial growth rate also needs to be elucidated. Finally, it will be fascinating to develop methods for studying the growth of bacterial colonies in 3D 38 , for which there are very few data at present 39 . Arriving at a predictive understanding of bacterial colony growth will constitute an advance not only in fundamental active matter physics, but also in practical fields such as food safety 40 .

Methods
Experimental methods. Overnight cultures of E. coli K-12 strain MG1655 were grown in M9 medium at 37°C and shaken at 200 rpm, harvested in exponential phase, and diluted 100-fold (to an optical density of ≈0.2-0.3 at 600 nm). A number of 1.5-2 μL droplets were then inoculated onto a thin layer of 2% M9 agarose filling a polymeric frame (Thermo Scientific Gene Frame AB-0630, ≈0.25 mm thick) stuck onto a standard microscope slide, and a cover slip was sealed on top. Growing colonies were observed at 37°C in phase contrast mode using a Nikon Eclipse E-800 microscope with a Nikon Plan Apo λ 100×/1.45 Ph3 Oil objective. Images were taken (Q-imaging Retiga 2000R camera) at regular intervals (typically 1 min). Using a published protocol 41 we obtained the timedependent centre of mass coordinates r i (t), length L i (t), and orientational unit vector ν i (t) of each cell i. The cell diameter, W ≈ 0.9 μm, varies little between cells or during growth. The length varies from~2 μm just after division to~4 μm just before. We report data from 32 colonies.
Image recording and processing. A schematic of the experimental set-up is shown in Supplementary Fig. 1a. For observations under phase contrast illumination intensity was minimised and exposure time varied from 20 to 100 ms. In a separate experiment, we performed confocal microscopy using fluorescent bacterial cells and agarose. For this case, the dye used for agarose staining was rhodamine B (0.02% w/v) (red channel); this was used to enhance contrast with the GFP-labelled cells (green channel) and the glass that is the black area in the images. Rhodamine B preferentially concentrates at the agarose surface. Sideways reconstruction of an image stack, Supplementary Fig. 1b, shows that cells are embedded just within the agarose beneath the glass side.
Analysis of phase contrast images was performed using ImageJ and Matlab. Initially, frames were cropped to the actual area covered by the colonies at the final stage (the field of view achievable is at least three times larger). The frames were then adjusted in brightness and contrast, and corrected for inevitable drifts (a few) in the agarose pad, especially toward the beginning of the imaging period. The latter was due to the agarose taking~10 min to cool to room temperature, having been poured into the frame at its melting point. The shift correction was performed using ImageJ (macro StackReg 42 ). Finally, the frames were segmented using the Matlab-based Schnitzcells software 41 to extract the positions, orientations and lengths of the individual cells. In order to extract velocities, the trajectories were stitched together in MATLAB using the standard Crocker and Grier code 43 .
Active field theory simulations. Numerical simulations of Eqs. (1), (11), and (21) were performed by using finite difference methods to discretise Eqs. (1) and (21). For reasons of numerical stability, we have set the growth rate Λ to 0 when the density was locally larger than 2ρ 0 . Note also that in Eq. (21) we need v, rather than ρv which is the natural variable in Eq. (11): to avoid problems arising from division by a very small number we have set v = 0 when ρ was below a threshold (ρ 0 /2 worked well for our parameter set).
In wet active gels, the key control parameter to determine when active turbulence sets in is aL 2 s K , where a is the active stress parameter (in our case equal to a 0 ρ), L s is the system size, and K is the liquid crystalline elasticity 6 . Active turbulence sets in when this control parameter exceeds a critical value: therefore an infinite system always develops active turbulence. On the contrary, for dry active fluids the turbulent regime sets in when a 0 ρ γK exceeds a threshold 35 . This holds for our growing system as well, and we find the threshold is~2, as can be shown in Supplementary Figs. 8, 9, where we plot the density of points with low orientational order. This plot also shows that the density of low order regions, which we assume should scale as l À2 a scales linearly with a 0 close to the transition to active turbulence, corresponding to a scaling of l a~( a 0 + C) −1/2 in this regime (this is consistent with the behaviour found in the ref. 35 ). More data are required to firmly prove this scaling, as the range of a 0 we investigated is relatively narrow.
The fact that in our experiments there are microdomains within our colonies suggests we are in the regime where nematic order is disrupted and defects arise (active turbulence). Therefore a suitable regime to describe experiments within our active field theory requires a 0 ρ γK ! 2.
Code availability. The Matlab codes used for the analysis and in-house code written for active field theory simulations are available from the corresponding author upon request.

Data availability
Data on bacterial colony growth are available from the corresponding author upon request. The tracked colony data used in this manuscript are also available at https://dx.