# Defect states in compressible active polar fluids with turnover

Ludovic Dumoulin, Carles Blanch-Mercader, Karsten Kruse

► To cite this version:

Ludovic Dumoulin, Carles Blanch-Mercader, Karsten Kruse. Defect states in compressible active polar fluids with turnover. 2025. hal-05393631

HAL Id: hal-05393631

https://hal.science/hal-05393631v1

Preprint submitted on 2 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Defect states in compressible active polar fluids with turnover

Ludovic Dumoulin,¹,² Carles Blanch-Mercader,³,\* and Karsten Kruse¹,²,†

¹Department of Biochemistry, University of Geneva, 1211 Geneva, Switzerland

²Department of Theoretical Physics, University of Geneva, 1211 Geneva, Switzerland

³Physico-Chimie Curie, Curie Institute, PSL University, 75005 Paris, France

(Dated: June 11, 2025)

Biological active matter like the cytoskeleton or tissues are characterized by their ability to transform chemical energy into mechanical stress. In addition, it often exhibits orientational order, which is essential for many cellular and morphogenetic processes. Experimental evidence suggests that defects in the orientational order field play an important role in organizing active stress. However, defects tend to annihilate unless the material is in a chaotic state or hydrodynamic interactions are suppressed. Using a hydrodynamic description of compressible active polar fluids, we show that turnover readily leads to a stabilization of defects. Depending on the turnover rate, topological defects arrange in a multitude of different phases, including lattices, active foams, and vortex glasses. Our work suggests that turnover plays a crucial role for organizing biological active matter.

## I. INTRODUCTION

The physics of active matter [1] is known to play a key role in various morphogenetic processes that occur in single cells or tissues. Well-known examples include the formation of lamellipodia in migrating single cells [2], apical constriction during gastrulation of the fruit fly *Drosophila melanogaster* [3], or cortical actin flows during head-tail symmetry breaking in the nematode *Caenorhabditis elegans* [4]. These processes are driven by gradients in active stress, which is generated autonomously by the biological systems through the conversion of chemical energy released during the hydrolysis of adenosine-triphosphate (ATP) into mechanical work. How gradients in active stress are organized to result in the observed morphogenetic events remains largely unknown.

On a molecular level, active stress is generated by motor proteins that interact with protein filaments, for example, myosin motors and actin filaments. The stress generated in this way is anisotropic due to the elongated form of the filaments. This molecular anisotropy can translate to macroscopic orientational order when filaments align on cellular scales or when elongated cells align on tissue scales. The corresponding orientational order fields are often nematic [5, 6], when the arrangement is invariant under inversion of the orientation, or polar [7–9], when this invariance is broken. Higher order orientational order is also possible [10]. Active stress is tightly coupled to the orientational order in a cell or tissue [11].

In recent years, experimental evidence has accumulated that topological defects play an important part in organizing active stress [12–15]. These defects are singularities of the orientational field and are characterized by their topological charge [16]. Defects can themselves be oriented [17, 18], but the interaction of defects is not fully

accounted for by local properties [19]. Topological defects can lead to the accumulation of stress [12, 20], which could be related to the localization of cell death [12], accumulation [13] as well as differentiation [15], but also to morphogenetic processes like the formation of protrusions [14, 15, 21, 22].

In passive materials with orientational order, topological defects tend to spontaneously annihilate, unless their existence is required by the boundary conditions or the system's topology through the Poincaré-Hopf theorem. In contrast, active stress can lead to a continuous proliferation of singularities in the orientation field [5, 23, 24]. However, this proliferation is related to chaotic dynamics – or low Reynolds number turbulence [25] – such that it does not seem to be appropriate for a robust organization of active stress. Mechanisms for stabilizing topological defects in active matter remain largely elusive. Theoretical studies suggest that anisotropic substrates [26] and coupling to extrinsic curvature can generate orientational defect alignment [27]. In turbulent states, defect binding has been reported [28]. Finally, substrate friction leads to stable defect lattices [29–31].

In spite of a large body of theoretical work on the dynamics of defects [32] and the accumulating evidence that topological defects organize active stress during developmental processes [33, 34], the effects of turnover on defects in active fluids are hardly known. In this work, we introduce generalized hydrodynamic equations for compressible active polar fluids in the presence of turnover, Fig. 1(a). We find that turnover can stabilize topological defects. Using long-time numerical solutions of the dynamic equations in large systems, we find a number of different states, including active foams, density waves, and vortex glass states, Fig. 1(b-d), Movies 1-3. In the presence of anisotropic active stress additional phases appear. We show that the linear dynamics of density and polarity perturbations rationalize selection of the phases. These results indicate a central role of turnover in stabilizing topological defects and suggest a basic mechanism for a well-controlled organization of active stress.

<footnote>\* carles.blanch-mercader@curie.fr</footnote>

<footnote>† karsten.kruse@unige.ch</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Compressible active polar fluids in the presence of turnover. (a) Schematic of the physical processes considered in this work: turnover, active stress, and density-polarity coupling. (b-d) Snapshots of the density heatmaps from numerical solutions of the dynamic equations presented in Sec. II. Parameters are $ρ_0 = 0.45$ (b), $ρ_0 = 0.6$ (c), and $ρ_0 = 0.75$ (d) as well as $τ = 5$ and $ζ_ρ = 4$. In (b-d), the system size is $L = 50$. Other parameter values as in Table I. Maximal and minimal densities are, respectively, $ρ_{max} = 0.96$ and $ρ_{min} = 0.37$ (b), $ρ_{max} = 0.87$ and $ρ_{min} = 0.44$ (c), $ρ_{max} = 0.79$ and $ρ_{min} = 0.58$ (d).</figcaption></figure>

## II. HYDRODYNAMICS OF AN ACTIVE POLAR FLUID IN THE PRESENCE OF TURNOVER

Motivated by the actin cortex, we study a two-component system consisting of a compressible active polar fluid and a passive fluid on a two-dimensional space [35–38], Fig. 1(a). In the context of the actin cortex, the former represents the actomyosin network and the latter the cytoplasm containing actin monomers and

unbound molecular motors. In the following, we derive hydrodynamic equations that describe the system’s behavior on large space and time scales. As in Ref. [37], we neglect permeation flows of the passive component and focus only on the active component from here on. We denote the corresponding mass density by $ρ$ and its velocity field by $v$. In addition, the active component can exhibit macroscopic polar order, which is captured by the polarity field $p$.
The equilibrium behaviour of the system is determined by the total free-energy $\mathcal{F} = \int f d\mathbf{r}$, where $d\mathbf{r}$ is an infinitesimal area element and $f$ is the free-energy density $f$ given by

$$f = \frac{a}{4}\rho^4 + \rho^2 \left( -\frac{\chi}{2} \frac{\rho}{\rho_0} |\mathbf{p}|^2 + \frac{\chi}{4} |\mathbf{p}|^4 + \frac{\kappa}{2} (\nabla \mathbf{p})^2 \right). \quad (1)$$

The first term captures the energetic cost for concentrating the active component with an elastic coefficient $a$. We chose this term to be quartic in the mass density $\rho$ to prevent density divergences in the active case. As long as it is large enough, the precise value of the exponent is not important, though. The remaining terms correspond to the Landau-Ginzburg free-energy for the polarity field $\mathbf{p}$ [39]. For $\chi > 0$, we have $|\mathbf{p}|^2 = \rho/\rho_0$ in equilibrium, reflecting an increase of polar order with increasing density. For the time being, $\rho_0$ is some reference density. The final term penalizes distortions of the polarity field with an elastic constant $\kappa$ [39]. We neglect a potential spontaneous splay term.

In the presence of polar order, the hydrostatic pressure generalizes to the Ericksen stress $\sigma^e$ with components

$$\sigma_{\alpha\beta}^{e} = (f - \rho\mu)\delta_{\alpha\beta} - \frac{\partial f}{\partial(\partial_{\beta}p_{\gamma})}\partial_{\alpha}p_{\gamma}, \quad (2)$$

where greek indices run from 1 to 2 and $\delta_{\alpha\beta} = 1$ if $\alpha = \beta$ and zero otherwise. We adopt the convention that summation is implied over repeated indices. In the above expression, $\mu = \frac{\delta\mathcal{F}}{\delta\rho}$ is the active fluid's chemical potential and $h_\alpha = -\frac{\delta\mathcal{F}}{\delta p_\alpha}$ are the components of the molecular field. We give their expressions in App. A. The isotropic part of $\sigma^e$ can be interpreted as a pressure, $\Pi^{\text{eff}} = -f + \rho\mu$ with

$$\Pi^{\text{eff}} = \frac{3}{4} a \rho^4 - \rho^2 \left( \chi \frac{\rho}{\rho_0} |\mathbf{p}|^2 - \frac{\chi}{4} |\mathbf{p}|^4 - \frac{\kappa}{2} (\nabla \mathbf{p})^2 \right). \quad (3)$$

In equilibrium, when $|\mathbf{p}|^2 = \rho/\rho_0$, the effective isothermal compressibility $\frac{1}{\rho}\frac{\partial\rho}{\partial\Pi}$ is positive only if $\rho_0 > \sqrt{5\chi/6a}$. For smaller values of $\rho_0$, the system is thermodynamically unstable.

We now turn to the dynamics of the system. We assume that it is maintained out of equilibrium by the continuous consumption of a chemical fuel. In the case of the actin cortex, this fuel is adenosine tri-phosphate (ATP), which is hydrolyzed to form adenosine di-phosphate (ADP) and inorganic phosphate ($P_i$). The hydrodynamic description is based on conservation laws for mass and momentum and on the broken continuous symmetry under rotations resulting from the polar field [40].

Mass conservation of the active fluid reads

$$\partial_t \rho + \partial_\alpha (\rho v_\alpha) + \partial_\alpha j_\alpha = -\tau^{-1}(\rho - \rho_0). \quad (4)$$

The total material current is decomposed into a convective part $\rho\mathbf{v}$ and a diffusive part $\mathbf{j}$. The source terms account for assembly and disassembly of the active fluid.

We restrict ourselves to small deviations from the target density $\rho_0$ and denote by $\tau$ the characteristic time scale of the disassembly process. The values of $\rho_0$ and $\tau$ are related to the fluid's assembly and disassembly rates $k_a$ and $k_d$ via $\tau = k_d^{-1}$ and $\rho_0 = k_a/k_d$. In the context of tissues, these rates are associated with cell division and cell death processes, [41, 42]. The densities of ATP, ADP, and $P_i$ are not explicitly included in this description, as they are assumed to be spatially and temporally constant.

As cytoskeletal flows or tissue flows are typically associated with a low Reynolds number, we neglect inertial effects. From now on, $\rho$ denotes the number density of the active fluid, which is related to its mass density by some molecular mass that we consider to be constant. Furthermore, momentum conservation reduces to the force balance equation

$$\partial_{\beta}\sigma_{\alpha\beta}^{\text{tot}} = \xi v_{\alpha} \quad (5)$$

where $\sigma^{\text{tot}}$ denotes the total stress tensor. The right hand side corresponds to a viscous friction term with friction coefficient $\xi$. For the actomyosin cortex, it mainly results from interactions between the actomyosin network and the cell membrane or an underlying substrate although permeation of the actomyosin network by the cytosol also contributes to this term.

To close the system of equations, we need to provide constitutive equations of the active fluid. Following the standard approach of nonequilibrium thermodynamics [43], we identify the thermodynamic forces conjugate to the fluxes and express the latter to linear order in terms of the former. Explicitly, the fluxes are the symmetric part of the deviatoric stress $\sigma^d \equiv \sigma^{\text{tot},s} - \sigma^{e,s}$, the diffusive current $\mathbf{j}$, the convected co-rotational derivative of the polar field $D\mathbf{p}/Dt$, and the rate $r$ of ATP hydrolysis. Here, the superscript $s$ denotes the symmetric part of the respective tensors and the convected co-rotational derivative is defined by $Dp_\alpha/Dt = \partial_t p_\alpha + v_\beta \partial_\beta p_\alpha + \omega_{\alpha\beta} p_\beta$, where $\omega_{\alpha\beta} = (\partial_\alpha v_\beta - \partial_\beta v_\alpha)/2$ are the components of the vorticity tensor. Note that the total stress tensor in Eq. (5) also contains an anti-symmetric part such that

$$\sigma_{\alpha\beta}^{\text{tot}} = \sigma_{\alpha\beta}^{d} + \sigma_{\alpha\beta}^{e} + \frac{1}{2}(p_{\alpha}h_{\beta} - p_{\beta}h_{\alpha}). \quad (6)$$

The conjugated thermodynamic forces are the symmetric part of the velocity gradient tensor with components $v_{\alpha\beta} = (\partial_\alpha v_\beta + \partial_\beta v_\alpha)/2$, the gradient of the active fluid's chemical potential $\mu$, the molecular field $\mathbf{h}$, and the difference in chemical potentials of ATP and its hydrolysis products $\Delta\mu = \mu_{\text{ATP}} - \mu_{\text{ADP}} - \mu_{\text{P}}$. Explicitly, the constitutive equations are

$$\sigma_{\alpha\beta}^{d} = 2\eta v_{\alpha\beta} + \frac{\nu}{2}(h_{\alpha}p_{\beta} + h_{\beta}p_{\alpha}) + \rho\Delta\mu (\zeta_{\rho}\rho^{2}\delta_{\alpha\beta} + \zeta_{p}p_{\alpha}p_{\beta}) \quad (7)$$

$$j_{\alpha} = -\gamma \partial_{\alpha} \mu \quad (8)$$

$$\frac{Dp_{\alpha}}{Dt} = -\nu p_{\beta} v_{\alpha\beta} + \frac{1}{\Gamma} h_{\alpha}. \quad (9)$$
Here, for simplicity, among others we have neglected a term proportional to **p** in Eq. (8) and terms containing gradients of the chemical potential $\mu$ in Eqs. (7) and (9).

The constitutive equations contain a number of phenomenological parameters that characterize the material properties. The shear viscosity $\eta$, the rotational viscosity $\Gamma$, the mobility coefficient $\gamma$, and the flow-alignment coefficient $\nu$ are familiar from conventional liquid crystals. In addition, there are active coefficients $\zeta_{\rho}$ and $\zeta_{p}$ that couple the deviatoric stress $\sigma^{\mathrm{d}}$ to the chemical potential difference $\Delta\mu$. We have neglected a potential coupling between $\Delta\mu$ and the time derivative of the polarization as it amounts to a redefinition of the equilibrium polarization amplitude [44]. In agreement with the Curie principle, we have used the polarization **p** to couple fluxes and forces of different tensorial order. Furthermore, we included a dependence on the density $\rho$ in the active terms, because active stress is expected to increase with the active fluid's density. Similar hypotheses have been used in other descriptions of the acto-myosin cortex [37, 45].

In the following, we use the dynamic equations in non-dimensionalized form. Unless stated otherwise we use the parameter values given in Table I. The values of $\tau$, $\rho_0$, and $\zeta_{\rho}$ are varied. Additionally, in Eq. (7) we ignored an isotropic viscous stress $\bar{\eta}\nu_{\gamma\gamma}\delta_{\alpha\beta}$, an isotropic active stress $\bar{\zeta}_{\rho}p_{\gamma}p_{\gamma}\delta_{\alpha\beta}$, and an isotropic elastic stress $\bar{\nu}h_{\gamma}p_{\gamma}\delta_{\alpha\beta}$, as we checked that our results are robust to variations in the parameters $\bar{\zeta}_{\rho}$ and $\bar{\nu}$. Finally, we set $\Delta\mu = 1$ to simplify the notation.

### III. DEFECT STABILIZATION THROUGH TURNOVER

Through minimization of the free energy associated with orientational order, pairs of topological defects with opposite charges tend to merge and annihilate upon relaxation to thermodynamic equilibrium [17–19]. In the absence of noise, this relaxation process stops, when the defect configuration with minimal energy and consistent with the boundary conditions is reached. In contrast, active processes can promote the separation of defect pairs and even increase their number [5, 46]. For active nematics, this requires for instance an anisotropic active stress [47]. In polar fluids, we find that turnover can stabilize topological defect pairs both in the presence or absence of active stress, Fig. 2(a). In this section, we analyze the mechanism leading to defect pair separation. To this end, we first consider an isolated defect and then turn to the case of two defects.

#### A. An isolated defect

Our numerical results indicate that the fluid density increases with the distance from a defect center, Fig. 1(b,d). To approximate this continuous variation, we separate space into two domains. Placing the defect

<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Number of defects as function of time for different values of $\tau$ and $\rho_0 = 1$, $\zeta_{\rho} = 4$, $L = 10$. Solid lines: number of defects averaged over 50 runs per $\tau$, shaded area: standard deviation. (b) Snapshot of density ( heatmap) and polarity (black arrows) maps for a pair of topological defects. Parameters: $\rho_0 = 0.6$, $\tau = 1$, $\zeta_{\rho} = 12$. The image was cropped from a square domain of length $L = 10$. (c) Critical distance below which a pair of oppositely-charged defects annihilates, as function of $\rho_0$ and $\zeta_{\rho}$ for $\tau = 1$ and $L = 10$. The black region corresponds to a critical distance larger than half of the system size $L/2$. (d) Snapshot of density ( heatmap) and velocity (white arrows) maps for the same defect pair shown in panel (b).</figcaption></figure>

center at the origin $r = 0$, we define an inner circular domain with a radial distance $r < R$, where we take $\mathbf{p} = \mathbf{p}_{\mathrm{in}}$ with $|\mathbf{p}_{\mathrm{in}}|^2 = 0$ and $\rho = \rho_{\mathrm{in}} = \mathrm{const}$. This inner domain represents the defect core and $R$ represents the radius of a defect core. In the outer domain with $r > R$, we take $|\mathbf{p}|^2 = |\mathbf{p}_{\mathrm{out}}|^2$ and $\rho = \rho_{\mathrm{out}} = \mathrm{const}$. Since the outer domain is much larger than the inner domain, we assume that the density there is set by turnover and choose $\rho_{\mathrm{out}} = \rho_0$. In addition, we consider that the polar order parameter attains its equilibrium value, such that $|\mathbf{p}_{\mathrm{out}}|^2 = 1$. Here, for simplicity, we have chosen to neglect the effect of the Frank free energy, i.e., $\kappa = 0$. As a consequence, the Ericksen stress (6) reduces to the pressure (3). *A priori*, the value of $R$ is determined by the steady state polarity profile. However, since we neglect the Frank free energy, $R$ remains undetermined and is treated as a parameter in the following. In addition, as a consequence of this approximation, we cannot enforce chemical and mechanical balance simultaneously.

If there is a difference between $\rho_{\mathrm{in}}$ and $\rho_{\mathrm{out}}$, there can be a current across the interface between the two domains. This current is driven by differences in the stress and the chemical potential across the domain boundary. Given that the mobility coefficient $\gamma = 10^{-4}$ is small in our computations, we neglect the difference in chemical potential and consider only the stress difference as a driver of the flow. Such a flow can be sustained in the presence of turnover.

The advective flux across the interface $r = R$ is given
<table><caption>TABLE I. Values of the dimensionless parameters used in this work. Dimensionless parameters are denoted by the same symbol as the corresponding dimensional parameter, with a prime. The relation between dimensionless and dimensional parameters, as well as the choice of unit system, is detailed in Appendix B. * $a' = 4/3$ if $\zeta_\rho = 0$. For simplicity, we omit the prime notation in Secs. III-VI and in the figures of the main text.</caption><thead><tr><th>parameter</th><th>ρ'<sub>0</sub></th><th>a'</th><th>χ'</th><th>κ'</th><th>ν'</th><th>ζ'<sub>ρ</sub></th><th>ζ'<sub>ρ</sub></th><th>γ'</th><th>Γ'</th><th>τ'</th><th>Δμ</th></tr></thead><tbody><tr><td>value</td><td>0.4-1.5</td><td>4/3ζ'<sub>ρ</sub>*</td><td>0.1</td><td>10<sup>-4</sup></td><td>0</td><td>0-24</td><td>0</td><td>10<sup>-4</sup></td><td>1</td><td>0.01-100</td><td>1</td></tr></tbody></table>

by $v(R^-)\rho_{in}$, where $v(R^-)$ represents the radial component of the fluid velocity in the inner domain evaluated near the interface. To estimate this velocity, we approximate the fluid flow as a Darcy flow, which is proportional to the total stress difference across the interface, $v(R^-) = \Delta\sigma/(\ell\xi)$ with $\Delta\sigma = \sigma(R^+) - \sigma(R^-)$, where $\sigma(R^\pm)$ denote the total stress in the inner (-) and outer domains (+) evaluated near the interface. The parameter $\ell$ is the thickness of the interface, which we assume to be proportional to $R$. Combining Eqs. (2), (3), (6) and (7), the total stress difference $\Delta\sigma$ is given explicitly by

$$ \Delta\sigma = \frac{3}{4}\chi\rho_0^2 - \frac{3}{4}a(\rho_0^4 - \rho_{in}^4) + \zeta_\rho(\rho_0^3 - \rho_{in}^3). \quad (10) $$

Note that mechanical equilibrium $\Delta\sigma = 0$ imposes a density difference $\rho_{in} \neq \rho_0$. To make further progress, we assume the density inside the defect to be $\rho_{in} = \rho_0(1 - \varepsilon)$, where $\varepsilon \ll 1$, and expand (10) up to linear order with respect to $\varepsilon$,

$$ \Delta\sigma = \frac{3}{4}\chi\rho_0^2 - 3(a\rho_0 - \zeta_\rho)\rho_0^3\varepsilon + \mathcal{O}(\varepsilon^2) \quad (11) $$

To determine the value of $\rho_{in}$, we impose mass conservation: the rate at which the active fluid is produced in the inner region must equal the rate at which it flows out through the interface. Explicitly, total mass balance condition is given by

$$ \tau^{-1}(\rho_0 - \rho_{in})\pi R^2 = 2\pi R v(r = R^-)\rho_{in}. \quad (12) $$

Up to first order in $\varepsilon$, we then get from Eq. (12) the fluid velocity

$$ v(R^{-}) = \frac{\tau^{-1}R}{2}\varepsilon, \quad (13) $$

where $\varepsilon$ takes the form

$$ \varepsilon = \frac{3\chi\rho_0^2}{3\chi\rho_0^2 + 12\rho_0^3(a\rho_0 - \zeta_\rho) + 2\xi\tau^{-1}R\ell}. \quad (14) $$

Equation (13) shows that radial flows can arise as an interplay between the density variations that are generated near topological defects and turnover.

Finally to check, whether the result above is consistent with our approximation to neglect the diffusive flux across the interface, we estimate the diffusive flux $j = -\gamma(\mu_{out} - \mu_{in})/\ell$ and find up to lowest order in $\varepsilon$

$$ j = \gamma\rho_0\chi/\ell. \quad (15) $$

For the advective current we get up to the same order

$$ v(R^{-})\rho_{in} = 3\chi\rho_{0}^{2}/(4\xi\ell), \quad (16) $$

such that for the parameters in Table I, we have the ratio between these two fluxes is

$$ |j/(v(R^-)\rho_{in})| = 4\xi\gamma/(3\rho_0) \approx 10^{-4} \ll 1. \quad (17) $$

## B. A defect pair

In the case of two defects, the fluid outflow generated near the core of one defect can advect the other and vice versa, which tends to separate the two defects. On the contrary, elastic interactions between the polar fields of the two defects tend to attract defects. The scale of these interaction forces in two dimensions is $\sim \kappa/d$, [16], where $d$ is the distance between the defect centers. Assuming an overdamped dynamics for the defect positions with an effective friction coefficient $\bar{\gamma}$, [47, 48], one gets by balancing the above effects

$$ \dot{d} = 2v(r = d) - \frac{\kappa}{d\bar{\gamma}}. \quad (18) $$

where $v(r = d)$ represents the velocity field generated by a single defect at position $r = d$.

For small $d \ll R$, elastic interactions dominate, causing the defects to eventually merge. For large $d \gg R$, elastic interactions dominate again, because they are long-ranged, whereas the radial flows are localized near the defect center. However, at intermediate distances $d \sim R$, radial flows can overcome the elastic attraction, leading to defect separation. This occurs when $\dot{d} > 0$ at $d \sim R$, which implies that $\tau^{-1}R^2\varepsilon > \kappa/\bar{\gamma}$, where $\varepsilon$ is given by Eq. (14) and we used Eq. (13) to estimate the fluid velocity.

This expectation is confirmed by numerical integration of the dynamics of two oppositely charged defects, which qualitative agrees with the analytic considerations, Fig. 2(b-d), Movies 4 and 5. Numerical results show that turnover of the active fluid can stabilize topological defects, Fig. 2(c). In large systems, this leads to a sustained finite density of defects, Fig. 3, maintained up to critical values of the density $\rho_0$ and the activity $\zeta_\rho$. This is compatible with our simplified analysis above: as $\rho_0$ and $\zeta_\rho$ increase, the value of $\varepsilon$ (14) decreases, leading to a reduction in the amplitude of the radial flows (13), and thus repulsive interaction between defects. Furthermore, the energetic cost of a defect increases with increasing
the target density $\rho_0$, promoting defect annihilation. Finally, the defect density decreases with $\tau$, Fig. 2(a) and Fig. 3, showing that sufficiently fast turnover is necessary for defect stabilization, which is again consistent with our simplified analysis.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Defect density in the $(\rho_0, \zeta_\rho)$-plane for different values of $\tau$. For each data point, a single numerical solution for a system size of $L = 10$ has been analyzed. White dots indicate the critical activity $\zeta_\rho^c$ obtained by solving Eq. (21) with $a = 4\zeta_\rho^c/3$ (values of $\zeta_\rho^c < 1$ are not shown). Blue symbols indicate solutions for $\rho_0 = 0.4$ (diamond), 0.8 (triangle), and 1.2 (circle) shown in Fig. 5 below. Other parameter values as in Table I.</figcaption></figure>

## IV. DEFECT STATES FOR ISOTROPIC ACTIVE STRESS

In this section, we show how topological defects arrange as a function of the target density $\rho_0$, which is set by the turnover rates. We start with a linear stability analysis of the homogenous polarized state. Then, we present numerical solutions of the dynamic equations, and finally rationalize the observed sequence of states using a linear stability analysis of the isotropic state.

### A. Linear stability analysis

Consider the uniform polarized state, $\rho = \rho_0$ and $p = p_0$ with $|p_0|^2 = 1$, which is a steady state of the dynamic equations (1)-(9). The linearized dynamics for a small perturbation of this state is given in App. C. The ensuing stability diagram shows that this state is unstable for a critical level of active stresses $\zeta_\rho > \zeta_\rho^c$, where the critical value $\zeta_\rho^c$ is a nonmonotonic function of the target density $\rho_0$, Fig. 4(a).

The instability is mainly driven by fluctuations in the density. To see this, we computed the growth rate associated with density fluctuations $\lambda_\rho$ in the special case

<figure><img src="image_5.png" /><figcaption>FIG. 4. Linear stability analysis of homogenous states. (a) Stability diagram of the homogenous polar stationary state as a function of $\rho_0$ and $\zeta_\rho$ for $\tau = 5$ and $a = 4/3\zeta_\rho^c$. Solid lines: density and polarity perturbations, diamonds: only density perturbations. (b) Fastest growth rate of density (red) and polarity (blue) perturbations of the homogenous isotropic state for $\zeta_\rho = 4$ and $\tau = 5$, other parameters as in Table I.</figcaption></figure>

$\chi = 0$, which takes the form

$$ \lambda_{\rho}(q) = - \frac{(q^2 - q_c^2)^2}{\tau q_c^4 (1 + 2q^2)} + \frac{3q^2 \rho_0^3 (\zeta_{\rho} - \zeta_{\rho}^c)}{1 + 2q^2}. \quad (19) $$

The variable $q_c$ denotes the wavenumber of the perturbation with the maximal growth rate $\lambda_\rho$ at the critical value $\zeta_\rho = \zeta_\rho^c$ and is given by

$$ q_c = \frac{1}{(2\gamma A(\chi = 0)\tau)^{1/4}}. \quad (20) $$

where $A = 3a\rho_0^2 - 5\chi/2$, see App. C 1.

The variable $\zeta_\rho^c$ is the critical level of the active stress at which the growth rate at $q = q_c$ vanishes, i.e. $\lambda_\rho(q_c, \zeta_\rho = \zeta_\rho^c) = 0$, and is given by

$$ \zeta_{\rho}^{c} = \frac{A(\chi = 0)}{3\rho_{0}} + \left( \sqrt{\frac{A(\chi = 0)\gamma}{3\rho_{0}^{3}}} + \sqrt{\frac{2}{3\rho_{0}^{3}\tau}} \right)^{2} \quad (21) $$

This critical activity threshold provides a good approximation for the value $\zeta_\rho^c$, when the value of $\chi = 0.1$, Fig. 4(a).

Let us point out that the linear dynamics are qualitatively different in absence of turnover. In that case, density perturbations with $q = 0$ are marginal, implying a type II instability according to the classification of Cross and Hohenberg [49]. In contrast, the instability in the presence of turnover is of type I.

### B. Nonlinear dynamics

To explore the system dynamics beyond linear stability, we numerically solved the dynamic equations (1)-(9) using a custom code written in Julia [50]. Details are given in App. D. The code can be found at [51]. Typically, we used square domains of size $L = 10$ with a discretization length $\Delta x = 10^{-2}$ and a total simulated time
of $10^5$. Selected simulations were run for a total simulated time of $2 \cdot 10^6$ or square domains of size $L = 50$. Units are set as described in App. B.

Simulations are initiated with a slight perturbation of the uniform disordered state by randomly perturbing both the uniform density, $\rho = \rho_0$, and the uniform polarity field, $\mathbf{p} = \mathbf{0}$. We observed a variety of different asymptotic behaviors. In the following, we describe the different phases observed in the computations, Fig. 5.

### 1. Puncta and networks

For small target densities $\rho_0$, $\rho_0 = 0.4$ in Fig. 5 and Movie 6, the system spontaneously forms puncta of elevated density within a background with a density around $\rho_0$. The majority of these puncta are connected to other puncta by lines of intermediate density that formed the edges of a network with puncta at its vertices. We denote the regions enclosed by edges as cells.

The network is dynamic and undergoes topological changes as puncta move, fuse, and split. Simultaneously, polar order is present for these parameter values, Fig. 5A. The polarity does not exhibit a preferred orientation at the cell edges. Furthermore, the polar field exhibits defects with integer charges $\pm 1$. Defect centers are mainly located within the cells.

### 2. Active foams

With increasing target density $\rho_0$, the density along the edges increases, $\rho_0 = 0.5$ in Fig. 5 and Movie 7. The system forms two sorts of cells: One kind contains a single integer topological defect of the polar order field, whereas cells of the other kind do not contain defects. In the latter, the polarity is typically aligned throughout the whole cell. Some edges have a vertex that is not connected to other edges. This typically occurs in cells without a defect in the polar order field. Similar to the case of $\rho_0 = 0.4$, this state is dynamic and we denote it in the following as an active foam. We discuss similarities and differences with active foam states reported in other systems in Sect. VII.

### 3. Density waves

As $\rho_0$ is increased further, traveling density waves appear, $\rho_0 = 0.6$ and $\rho_0 = 0.65$ in Fig. 5 and Movies 2 and 8. At the same time, uniform polar order emerges throughout the whole system. In contrast to lower densities, the cell shapes are elongated in one direction. The direction of elongation is on average perpendicular to the direction of motion of the density waves, but not correlated with the direction of the polarity field. These states are identical to the waves obtained when starting with a fully polarized state $\mathbf{p} = \mathbf{p}_0 = \text{const}$ with $|\mathbf{p}_0|^2 = 1$. At $\rho_0 = 0.7$,

Fig. 5 and Movie 9, we observe re-entrant behavior of active foams. This time the interfaces are thicker along the edges compared to active foams at lower densities, Fig. 5B. This state also features traveling density waves.

### 4. Vortex glass

For a target density of $\rho_0 = 0.75$, topological defects with cores of two different sizes emerge, Fig. 5 and Movie 3. The state is dynamic with “small” defects flowing between “large” ones. Due to the similarity of this state with states observed in the complex Ginzburg-Landau equation, we denote this state as a vortex glass state [52]. Note, however, that the disclination points with charge $+1$ can also be asters and spirals in addition to vortices, Fig. 5C.

### 5. Defect lattices

As the target density is increased further, the size of defect cores is again monodisperse. For densities around $\rho_0 = 0.8$, locally, regular defect lattices appear spontaneously, Fig. 5D and Movie 10. The topological defects have charges of $\pm 1$ and the lattices are either square or hexagonal. We discuss the lattice phase in more detail in Sect. V below.

### 6. Homogenous uniformly polarized state

For target densities $\rho_0$ beyond a critical value, the system asymptotically settles into the homogenous state with uniform polarity, Movie 11. Indeed, the critical value $\zeta_{\rho}^c$, Eq. (21), increases with the density $\rho_0$. By increasing $\rho_0$ and keeping all other parameter values constant, the homogenous uniformly polarized state becomes linearly stable. Our numerical results suggest that for sufficiently large densities $\rho_0$ the state is globally stable except for a uniform rotation of the polar field.

## C. Competition between density and polarity perturbations

We can rationalize the transitions between the states in Fig. (5) as a function of the average density $\rho_0$ by comparing the rates associated with the fastest growing mode of density fluctuations $\lambda_{\rho}^m$ and polarity fluctuations $\lambda_{\rho}^p$ around the isotropic homogeneous state. The calculation of these two rates is given in App. C. For the parameters used in Fig. 5, we present their values as a function of $\rho_0$ in Fig. 4(b).

For intermediate target densities $\rho_0 \gtrsim 0.35$ and $\rho_0 \lesssim 0.7$, density fluctuations grow faster than polarity fluctuations, Fig. 4(b). Correspondingly, density patterns dominate in this parameter regime and impose polarity
<figure><img src="image_6.png" /><figcaption>FIG. 5. Snapshots of the density ( heatmaps) and polarity ( black arrows) of asymptotic states of the dynamic equations (1)-(9) for different target densities ρ₀ for τ = 5 and ζρ = 4. The system size is L = 10. Other parameters as in Table I. The same states are presented for a system of length L = 50 in Fig. 10. The bottom right panel shows zoomed-in views of selected snapshots, labeled A–D. Maximal and minimal densities are indicated in legends.</figcaption></figure>

patterns, Fig. 5. With increasing average density ρ₀, the extent of the high-density domains increases.

Because the growth rate $\lambda_p^m$ grows quadratically with $\rho_0$, whereas the growth rate $\lambda_p^m$ is non-monotonic with $\rho_0$, there is a cross-over at about $\rho_0 \approx 0.75$ to a regime where polarity fluctuations grow faster than density fluctuations, Fig. 4(b). In this case, the final pattern results from the interaction between defects, see states at $\rho_0 = 0.75$ and $\rho_0 = 0.8$ in Fig. 5. At even larger tar-

get densities $\rho_0 > 0.85$, the system always relaxes to a polarized homogeneous phase, Fig. 5.

V. DEFECT LATTICES

In the defect lattices we described in the previous section, the topological defects in the polarity field locally arrange into regular polygons. This behavior is robust
<figure><img src="image_7.png" /><figcaption>FIG. 6. Examples of square and hexagonal defect lattices. (a,c) Density heatmaps at $t = 2 \cdot 10^6$ for $\rho_0 = 0.7$, $\tau = 1$, and $\zeta_\rho = 10$ (a) and for $\rho_0 = 1.3$, $\tau = 0.2$, and $\zeta_\rho = 1$ (c). (b,d) Angle of the polarity field corresponding to (a) and (c). System size $L = 10$, other parameters as in Table I. The corresponding time evolution of the shape order $\Gamma_p$ in Eq. (23) is given in Fig. 11.</figcaption></figure>

across varying parameter values. To further characterize these states, we perfromed longer simulations with a total time of $t = 2 \cdot 10^6$. Over time, a global crystalline order of defect positions emerged, Fig. 6. We observed both, square and hexagonal lattices, Fig. 6(a,c).

Let us recall that numerical computations start from a slightly perturbed uniform disordered state. Due to local spontaneous symmetry breaking, defects in the polarity field appear throughout the system, initially forming an irregular arrangement. The defects created in this way are stabilized by the mechanism discussed in Sect. III and eventually arrange in a nearly disclination-free lattice as in Fig. 6(a,c).

A striking feature of these lattices is that they move at constant velocity without defect rearrangements, Movies 12 and 13. Defect lattices move at constant velocity along a direction that is not determined by the symmetry axes.

A detailed inspection reveals an alternating pattern of defects with topological charges +1 and −1, Figs. 6(b,d) and 11(a,d). In the one-constant approximation, isolated +1 defects of either type - aster, spiral, and vortex - have the same total free-energy, [16]. These subtypes appear in both square and hexagonal defect lattices, but do not show any specific spatial organization.

The formation of defect lattices with clearly peaked

<figure><img src="image_8.png" /><figcaption>FIG. 7. Transition from square to hexagonal defect lattices. (a) Average of the shape function $\langle \Gamma_p - \Gamma_p^0 \rangle$ given by Eq. (23). For each value of $\rho_0$, the average is taken with respect to numerical solutions with $\tau = 0.1, 0.2, 0.5, 1, 1.2, 1.7, 2.5, 5, 10, 100$ and $\zeta_\rho = 0, 1, 2, 4, 6, \dots, 24$. Only solutions with defect density larger than 1.5 were analyzed. (b) Shape function $\Gamma_p - \Gamma_p^0$ as a function of the defect density for the solutions used in (a). Error bars in panel (a) are SD. In both panels, the domain size $L = 10$.</figcaption></figure>

structure factors takes too much time as if a thorough characterization of the corresponding phase space were feasible. To quantify the emerging crystalline order from shorter computer runs, we use the fact that the emergence of long-range translational order in the defect arrangement is preceded by long-range orientational order and use a method described in Ref. [10] for its quantification. We start by computing the Voronoi tesselation of defect centers, which assigns a polygonal cell to each defect based on its surrounding defects. From these polygons, we compute the shape function $\gamma_p$ for $p = 1, 2, 3, \dots$, which quantifies the degree of orientational order in the spatial arrangement of defect centers. Explicitly, for a single cell

$$ \gamma_p = \frac{1}{\Delta_p} \sum_{v=1}^{V} |\mathbf{r}_v|^p e^{ip\phi_v}, \qquad (22) $$

where $V$ is the number vertices of the cell, $\mathbf{r}_v$ the position of vertex $v$ relative to the defect center, $\phi_v$ its angle with the x-axis, and $\Delta_p = \sum_{v=1}^V |\mathbf{r}_v|^p$. For a regular $p$-polygon, one has $\gamma_p = 1$, because $|\mathbf{r}_v| = \text{const}$ and because the angle between two adjacent vertices, $\angle(\mathbf{r}_v, \mathbf{r}_{v+1})$, is $2\pi/p$. Finally, for each $p$, we average the absolute value of the shape function $\gamma_p$ over all defect centers,

$$ \Gamma_p = \langle |\gamma_p| \rangle. \qquad (23) $$

To assess the degree of order, we compare $\Gamma_p$ to $\Gamma_p^0$, its value for randomly distributed vertices, by computing the difference $\Gamma_p - \Gamma_p^0$.

The trajectories of $\Gamma_p - \Gamma_p^0$, which we computed for $p = 2, \dots, 6$, indicate early on, $t \lesssim 10^5$, the developing crystalline order present in the arrangement of defects, Fig. 11(c,f). For each parameter set that displays a defect
lattice, we use $\Gamma_p$ at $t = 1.4 \cdot 10^5$ to characterize the organization of the defects. Although, in principle, this quantity can depend on the initial conditions, we find that it essentially depends on the defect density Fig. 7b. Specifically, at sufficiently small defect densities, the dominant type of order is hexatic ($p = 6$). In contrast, at intermediate defect densities, the dominant order changes to quartic symmetry ($p = 4$), indicating that the emergence of quartic crystalline order is not an effect of packing constraints. This observation is further supported by defects with smaller cores arranging in triangular lattices, while those with bigger cores form square lattices. At intermediate defect densities ($> 1.5$ per unit area), the crystalline order in the spatial arrangement of defect centers can also be controlled by the target density $\rho_0$, Fig. 7a.

## VI. EFFECTS OF DYNAMIC ANISOTROPIES

In the previous section, we studied the effect of turnover in the absence of anisotropic active stresses and flow alignment. In the following, we show that flow alignment affects the size of defect cores and can induce a transition to the vortex glass state. Furthermore, we show the emergence of labyrinthine patterns in the presence of an anisotropic active stress.

### A. Flow alignment

The flow-alignment parameter $\nu$ controls the coupling between the dynamics of the polarity field and shear flows. In addition to affecting the magnitude of polarity, the value of $\nu$ also determines how shear flows reorient the polarity field.

To gain insight into the role of $\nu$, we first review a simpler case: a polarization field subjected to a linear shear flow of the form $\mathbf{v} = \gamma(-x, y)$, where $x$ and $y$ are cartesian coordinates and $\gamma > 0$ denotes the shear strength [16]. In this setting, compression occurs along the $x$-direction, while expansion occurs along the $y$-direction. For $\nu < -1$, which corresponds to rod-shaped conventional liquid crystals, the polarity tends to align with the expansive direction. In contrast, for $\nu > 1$, which is characteristic of disc-shaped liquid crystals, the polarity tends to align with the compressive direction.

In our computations, we observed that single topological defects can create an outward radial velocity field that decays with the distance from the defect center, as discussed in Section III. Consequently, compression occurs in the radial direction. Based on the reasoning above, we expect flow-alignment to favor vortex-type defects for $\nu < -1$, and aster-type defects for $\nu > 1$.

Next, we studied the effect of the flow-alignment parameter on defect lattices. Previously, for $\nu = 0$, we observed an alternating pattern of +1 and -1 defects, where the +1 defects could be of any type, aster, spiral, or vortex, Fig. 11. For $\nu \neq 0$, in contrast, flow alignment

leads to a preference of asters for $\nu > 0$ and of vortices for $\nu < 0$, Fig. 8. This is in agreement with our simplified analysis above. To assure continuity of the polarization field, the direction of the polarization field is opposite in neighboring asters or neighboring vortices, see arrows in insets in Fig. 8.

<figure><img src="image_9.png" /><figcaption>FIG. 8. Effect of flow-alignment on +1 topological defects. Density (heatmap) and polarity (black arrows) for $\rho_0 = 0.8$, $\tau = 0.2$, and $\zeta_\rho = 4$. Changing the flow alignment parameter ($\nu$) can turn all the +1 topological defects into vortices ($\nu = -1.1$) or asters ($\nu = 1.1$). In both cases, the fluid has not yet reached a steady state. The domain size is $L = 10$.</figcaption></figure>

Furthermore, changing the flow alignment parameter affects the density and distribution of defects. This is most readily observed at low target densities $\rho_0$, Figs. 8 and 12. The defect density is larger for $\nu < 0$ than for $\nu > 0$, Fig. 12. Also, it increases with $\nu$. At these intermediate stages of the solutions, the defects form strings.
<figure><img src="image_10.png" /><figcaption>FIG. 9. Effects of active anisotropic active stress. (a) Contractile anisotropic stress without flow alignment. Parameters: $ρ_0 = 0.5$, $τ = 4.0$, $ζ_ρ = 4$, $ν = 0$ and $ζ_p = -1$. (b) Extensile anisotropic stress with negative flow alignment. Parameters: $ρ_0 = 0.4$, $τ = 4$, $ζ_ρ = 4$, $ν = -0.5$ and $ζ_p = 1$. (c) Extensile anisotropic stress with positive flow alignment. Parameters: $ρ_0 = 0.6$, $τ = 1$, $ζ_ρ = 8$, $ν = 0.5$ and $ζ_p = 0.5$. (d) Extensile anisotropic stress without flow alignment. Parameters: $ρ_0 = 0.4$, $τ = 2$, $ζ_ρ = 8$, $ν = 0$ and $ζ_p = 1$. The domain size is $L = 10$. Other parameter values as in Table I.</figcaption></figure>

Furthermore, the density around a defect is anisotropic, see for example Fig. 12 case $ν = 0.5$ and $ρ_0 = 0.75$.

## B. Anisotropic active stress

Next we explore the effects of an anisotropic active stress $ζ_p \neq 0$ on the long-term system dynamics. Anisotropic active stress is known to qualitatively change the dynamics of an active polar fluid. For example, they lead to sustained spontaneous flows [40, 53–56].

In Figure 9(a) and Movie 14, we present a snapshot of the system for $ζ_p = -1$. We observe the emergence of a labyrinthine density pattern. The paths of the labyrinth can be closed, for example circular, or exhibit bifurcations. In regions free of bifurcations, the paths tend to align. Labyrinthine patterns have also been found in defect free active nematic fluids [57]. In contrast, to the patterns reported here, the labyrinths are formed by nematic domain walls, whereas the density is constant.

In Figure 9(b-d) and Movies 15-17, we present snapshots of solutions for extensile anisotropic active stress. In contrast to Fig. 9(a), we also consider flow alignment. For negative flow alignment, $ν < 0$, the fluid forms high density rings, some of which are connected by bridges, Fig. 9(b). Typically, the density is higher along segments that are not adjacent to segments of other rings.

For positive flow alignment, $ν > 0$, there are again large regions of low density, but they are limited by broader high-density regions, Fig. 9(c). Besides, there are several lines of high density between the low-density regions. In the absence of flow alignment, $ν = 0$, but still extensile anisotropic active stress, Fig. 9(d), we observe again an active foam state. Some adjacent foam bubbles have separate surfaces. All the structures shown in Fig. 9 are dynamic.

These examples of solutions in the presence of flow alignment and anisotropic active stress component show that they can either have rather subtle effects by affecting the type of defect or more dramatic consequences by leading to features not observed in their absence, namely, disordered high-density lines.

## VII. DISCUSSION

In this manuscript, we studied the effect of turnover on the dynamics of active polar fluids with a focus on couplings between the density and polarity fields. Through a combination of extensive numerical computations and analytical calculations, we found various macroscopic phases, most notably active foams and defect lattices. The former are characterized by low-density regions surrounded by high-density interfaces, forming a dynamic network that continuously re-arranges. The latter corresponds to arrays of integer topological defects that over time develop long-ranged crystal order, in particular triangular and square lattices. The topological defects are stabilized by an interplay between density gradients at the defect core and turnover. This mechanism is independent of the strength of the anisotropic active stress, unlike the mechanism for defect pair separation in active nematic fluids, [47, 48]. Finally, we discussed two extensions of the theoretical description to account for alignment by shear flows and anisotropic active stress.

Theoretically, a contractile network with some similarity to the active foam discussed above has been obtained in Ref. [58]. In that work, turnover of an active fluid is assumed to be regulated by a biochemical reaction network. However, in that description, the active fluid is isotropic and the contractile network might approach a stationary state. Finally, an active foam state has been reported for so-called Malthusian flocks in the limit of constant density [59]. There, the "bubbles" of the foam are delimited by line defects in the orientational order field.

Experimentally, an active foam state has been reported in three-dimensional microtubule-motor mixtures [60]. Motor accumulation at microtubule ends is an important feature of this state. As our theory does not distinguish between filaments and motors, it is not obvious how to relate it to the experiments on motor-microtubule networks. Furthermore, in the experiment, the accumulation of motors at microtubule ends rather leads to high motor-filament densities in asters, which is opposite to
the depletion of active fluid generated by the free energy density (1). We conclude that the mechanisms underlying active foam formation are different in the experiment [60] and our theory. More recently, cytoplasmic compartmentalization by microtubule asters has been reported in zebrafish embryos and frog extracts [61]. Phenomenologically, these states resemble the active foam state shown in Fig. 5 for $\rho = 0.7$. Whether our description appropriately captures the mechanism underlying compartmentalization in these systems remains to be seen.

There have been a number of theoretical works on defect lattices in active fluids [29–31, 62–65]. In two-dimensional nematic active fluids, defect lattices were notably reported in the case of high surface friction compared to internal friction [30, 31, 65]. In these works the active anisotropic stress was extensile. As we have shown in Sect. III, in our polar fluid, turnover leads to a repulsive interaction between defects in presence of contractile activity. At intermediate values of the target density, this repulsion can lead to the formation of defect lattices.

Turnover is a general feature of living matter, on subcellular as well as tissue scales. We thus expect our findings to be relevant in various biological situations. It should be noted, though, that in the absence of fluctuations, as it was considered in this work, topological defects are not created by the dynamics of the contractile active fluid. Instead, they are generated through the initial condition. To robustly yield a certain density of defects, turnover could be regulated biochemically. In future work, it will be interesting to study the effects of couplings between biochemical regulation and the dynamics of active fluids.

## Appendix A: Expressions for the chemical potential, the molecular field and the Ericksen stress

In this section, we give the explicit expressions of the chemical potential $\mu = \frac{\delta F}{\delta \rho}$, the molecular field $\mathbf{h} = -\frac{\delta F}{\delta \mathbf{p}}$, and the Ericksen stress $\sigma_{\alpha\beta}^e$ given by Eq. (2). The total free-energy is $F = \int f d\mathbf{r}$, where $d\mathbf{r}$ is an infinitesimal area element and $f$ is the free-energy density given by Eq. (1). We obtain

$$ \mu = a\rho^3 + \rho \left( -\frac{3}{2}\chi \frac{\rho}{\rho_0} p_{\gamma}^2 + \frac{1}{2}\chi(p_{\gamma}^2)^2 + \kappa(\partial_{\gamma}p_{\delta})^2 \right) \quad (A1) $$

$$ h_{\alpha} = \chi \rho^2 \left( \frac{\rho}{\rho_0} - p_{\gamma}^2 \right) p_{\alpha} + \kappa \partial_{\beta} (\rho^2 \partial_{\beta} p_{\alpha}) \quad (A2) $$

$$ \sigma_{\alpha\beta}^e = - \left[ \frac{3}{4} a \rho^4 - \rho^2 \left( \chi \frac{\rho}{\rho_0} p_{\gamma}^2 - \frac{1}{4} \chi (p_{\gamma}^2)^2 \right) \right] \delta_{\alpha\beta} \\ - \rho^2 \left[ \frac{1}{2} \kappa (\partial_{\gamma} p_{\delta})^2 \right] \delta_{\alpha\beta} - \kappa \rho^2 (\partial_{\alpha} p_{\gamma}) (\partial_{\beta} p_{\gamma}). \quad (A3) $$

## Appendix B: Dimensionless equations

Let us begin by defining the units used in our theoretical description. The time unit is set by an arbitrary time scale $T$, the length scale by $\lambda = \sqrt{\eta\xi^{-1}}$ and the viscosity scale by $\eta$. Combining these expressions, we define a density scale as $\tilde{\rho} = \eta T \lambda^{-2}$.

In this unit system, the dimensionless parameters are: $a' = a\tilde{\rho}^4 T \eta^{-1}$, $\chi' = \chi\tilde{\rho}^2 T \eta^{-1}$, $\kappa' = \kappa\tilde{\rho}^2 T \eta^{-1} \lambda^{-2}$, $\rho_0' = \rho_0\tilde{\rho}^{-1}$, $\tau' = \tau T^{-1}$, $\nu' = \nu$, $\zeta_{\rho}' = \zeta_{\rho}\tilde{\rho}^3 T \eta^{-1}$, $\zeta_{p}' = \zeta_{p}\tilde{\rho} T \eta^{-1}$, $\gamma' = \gamma\eta\lambda^{-2}\tilde{\rho}^{-2}$, and $\Gamma' = \Gamma\eta^{-1}$.

Finally, we set $\Delta\mu = 1$ to simplify the notation.

## Appendix C: Linear stability analysis of the homogeneous state

In this section, we present the linear stability analysis of two homogeneous states: the polarized state and the isotropic state. We focus on the special case with no flow alignment, $\nu = 0$ and no anisotropic stress, $\zeta_p = 0$.

### 1. Polarized homogeneous state

To study the linear stability of the polarized state, we consider without limiting generality that the polar field is aligned along the $x$-axis, $\mathbf{p}_0 \cdot \mathbf{e}_y = 0$. We set $\rho = \rho_0 + \delta\rho$, $\mathbf{p} = \mathbf{p}_0 + \delta\mathbf{p}$ with $|\mathbf{p}_0| = 1$, and $\mathbf{v} = \delta\mathbf{v}$, and derive the dynamic equations to first order in the perturbations $\delta\rho$, $\delta\mathbf{p}$, and $\delta\mathbf{v}$. By setting $\delta\rho = \hat{\rho}(q_x, q_y)e^{i(q_x x + q_y y)}$ and similarly defining the Fourier components $\hat{v}_x$ and $\hat{v}_y$ for the velocity perturbation $\delta\mathbf{v}$ and $\hat{p}_x$ and $\hat{p}_y$ for the polarity perturbation, we use force balance (5) to express $\hat{v}_x$ and $\hat{v}_y$ in terms of $\hat{\rho}$, $\hat{p}_x$ and $\hat{p}_y$. Explicitly,

$$ \hat{v}_x = \frac{-i q_x \rho_0}{2 q^2 + 1} \left[ (A - 3 \rho_0 \zeta_{\rho}) \hat{\rho} - \chi \rho_0 \hat{p}_x \right] - \frac{1}{2} \frac{i q_y q^2 \rho_0^2}{q^2 + 1} \kappa \hat{p}_y \quad (C1) $$

$$ \hat{v}_y = \frac{-i q_y \rho_0}{2 q^2 + 1} \left[ (A - 3 \rho_0 \zeta_{\rho}) \hat{\rho} - \chi \rho_0 \hat{p}_x \right] + \frac{1}{2} \frac{i q_x q^2 \rho_0^2}{q^2 + 1} \kappa \hat{p}_y. \quad (C2) $$

Here, as in the main text, $A = 3\rho_0^2 - \frac{5}{2}\chi > 0$ and $q^2 = q_x^2 + q_y^2$ denotes the squared magnitude of the wavevector $\mathbf{q} = (q_x, q_y)$.

Substituting the expressions for $\hat{v}_x$ and $\hat{v}_y$ into the linearized equations (1)-(9), we obtain, to linear order in perturbations, the following dynamical equations for the density and the polarization

$$ \partial_t \hat{\rho} = -[\tau^{-1} + B(q)] \hat{\rho} + C(q) \hat{p}_x \quad (C3) $$

$$ \partial_t \hat{p}_x = \frac{\chi \rho_0}{\Gamma} \hat{\rho} - \rho_0 \frac{2\chi + \kappa q^2}{\Gamma} \hat{p}_x \quad (C4) $$

$$ \partial_t \hat{p}_y = -\kappa \rho_0^2 q^2 \left( \frac{1}{\Gamma} + \frac{1}{4} \frac{q^2}{q^2 + 1} \right) \hat{p}_y \quad (C5) $$
with $B(q) = q^2 [\rho_0^2(A - 3\rho_0\zeta_\rho)/(1 + 2q^2) + \gamma A]$ and $C(q) = q^2\chi\rho_0 [\rho_0^2/(1 + 2q^2) + \gamma]$. Note that perturbations in the polarization component perpendicular to the x-axis $p_y$ always decay, as shown in Eq. (C5). The linear dynamics of perturbations in $\rho$ and $p_x$ are coupled in general, and the resulting expressions for the two associated growth rates are rather involved, such that we do not present them here. However, valuable insights can be gained by considering the special case $\chi = 0$, in which the linear dynamics of $\rho$ and $p_x$ decouple. In this case, perturbations in $p_x$ always decay, see Eq. (C4), while perturbations in $\rho$ evolve with a growth rate $\lambda_\rho$, which can be recast in the form

$$ \lambda_{\rho}(q) = - \frac{(q^2 - q_c^2)^2}{\tau q_c^4 (1 + 2q^2)} + \frac{3q^2 \rho_0^3 (\zeta_{\rho} - \zeta_{\rho}^c)}{1 + 2q^2}. \quad (C6) $$

The variable $q_c$ denotes the wavenumber of the perturbation with the maximal growth rate $\lambda_\rho$ at the critical value $\zeta_\rho = \zeta_\rho^c$ and is given by

$$ q_c = \frac{1}{(2\gamma A(\chi = 0)\tau)^{1/4}}. \quad (C7) $$

The variable $\zeta_\rho^c$ is the critical level of the active stress at which the growth rate at $q = q_c$ vanishes, i.e. $\lambda_\rho(q_c, \zeta_\rho = \zeta_\rho^c) = 0$, and is given by

$$ \zeta_{\rho}^{c} = \frac{A(\chi = 0)}{3\rho_{0}} + \left( \sqrt{\frac{A(\chi = 0)\gamma}{3\rho_{0}^{3}}} + \sqrt{\frac{2}{3\rho_{0}^{3}\tau}} \right)^{2} \quad (C8) $$

## 2. Isotropic homogeneous state

In our numerical integration of the equations (1)-(9), we initialize the system in the isotropic state with a small random perturbation. In Sect. IV C, to interpret the states observed in numerics, we use a linear stability analysis of the isotropic state. Here, we provide a detailed explanation of this analysis.

We set $\rho = \rho_0 + \delta\rho$, $\mathbf{p} = \delta\mathbf{p}$, and $\mathbf{v} = \delta\mathbf{v}$, and derive the dynamic equations to first order in the perturbations $\delta\rho$, $\delta\mathbf{p}$, and $\delta\mathbf{v}$. In Fourier space (see previous section), we use force balance (5) to express $\hat{v}_x$ and $\hat{v}_y$ in terms of $\hat{\rho}$, $\hat{p}_x$ and $\hat{p}_y$. Explicitly,

$$ \hat{v}_x = 3E \frac{iq_x}{2q^2 + 1} \hat{\rho} \quad (C9) $$

$$ \hat{v}_y = 3E \frac{iq_y}{2q^2 + 1} \hat{\rho} \quad (C10) $$

with $q^2 = q_x^2 + q_y^2$ and $E = (\zeta_\rho - a\rho_0)\rho_0^2$. Substituting the expressions above for $\hat{v}_x$ and $\hat{v}_y$ into the linearized equations (1)-(9), we obtain the following dynamical equations for the density and the polarization

$$ \frac{d}{dt} \hat{\rho} = - \left\{ \frac{1}{\tau} + 3q^2 \left[ a\rho_0^2\gamma - \rho_0 \frac{E}{2q^2 + 1} \right] \right\} \hat{\rho}. \quad (C11) $$

$$ \frac{d}{dt} \hat{\mathbf{p}} = \frac{\rho_0^2}{\Gamma} (\chi - \kappa q^2) \hat{\mathbf{p}}. \quad (C12) $$

As long as $\chi > 0$, the isotropic state is unstable to perturbations in the polarity field, with a growth rate increasing quadratically with the target density $\rho_0$.

In Fig. 4, we represent the rates associated with the fastest growing mode of density fluctuations $\lambda_\rho^m$ and polarity fluctuations $\lambda_\rho^p$ around the isotropic homogeneous state. These rates were calculated using (C11) and (C12). For example, the rate $\lambda_\rho^m = \chi\rho_0^2/\Gamma$ is quadratic in the target density $\rho_0^2$.

## Appendix D: Numerical methods

In this appendix, we provide some details on the numerical used to integrate (1)-(9). The code is available at [51].

Spatial derivatives were calculated by first transforming the dynamic fields to the reciprocal space using the two-dimensional Fast Fourier Transform (FFT) technique, then multiplication with the corresponding wave vectors and, finally, transforming back to real space using FFT. Temporal integration was executed through the explicit Euler forward method. The computational grid encompassed at least 1008x1008 discrete points, with a spatial discretization of $\Delta x = 10^{-2}$. The maximal temporal discretization is set to $\Delta t_{max} = 10^{-2}$ and $\Delta t$ is chosen at every time unit $\Delta t$ according to the following rule $\Delta t = \min\{\Delta t_{max}, 5 \cdot 10^{-2} \Delta x/\max(||v||), 1.25\Delta t_{prev}\}$ with $\Delta t_{prev}$ being the previous time step. These computations were performed on GPUs.

At each time step, the chemical potential $\mu$ Eq. (A1), the molecular field $h_\alpha$ Eq. (A2), and the total stress minus the viscous stress $\sigma_{\alpha\beta}^{\text{nv}} = \sigma_{\alpha\beta}^{\text{tot}} - 2\eta v_{\alpha\beta}$, were computed in real space. Then, the velocity field was computed by solving the force balance Eq. (5) in Fourier space. Finally, the temporal density and polarity fields were updated by time integration of Eqs. (4) and (9).

## Appendix E: Supplementary figures

In this Section, we provide additional figures.
<figure><img src="image_11.png" /><figcaption>FIG. 10. Snapshots of the density (heatmaps) for different target densities $ρ_0$ and for $L = 50$, $τ = 5$ and $ζ_ρ = 4$. Other parameters as in Table I. The cases of $ρ_0 = 0.45$ and $ρ_0 = 0.55$ are not represented in Fig. 5.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 11. Example of a square and hexagonal defect lattice. (a,d) Snapshot at $t = 2*10^6$ of the density (heatmaps) and polarity (black arrows). (b,e) Voronoi tessellation for the patterns in (a) and (d). (c,f) Time evolution of the shape order parameter $\Gamma_p$ in Eq. (23) for $p = 2$ to 6 for the cases in (a) and (d). In panels (a-c) $\rho_0 = 0.7$, $\tau = 1$, and $\zeta_\rho = 10$ and in panels (d-f) $\rho_0 = 1.3$, $\tau = 0.2$, and $\zeta_\rho = 1$. The system size was set to $L = 10$. Other parameters as in Table I.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 12. Snapshots of the density (heatmaps) of solutions of the dynamic equations (1)-(9) for different target densities $ρ_0$ and flow alignment parameter $ν$ and for $τ = 5$ and $ζ_ρ = 4$. The system size is $L = 10$. Other parameters as in Table I.</figcaption></figure>
Appendix F: Movies

In this section, we provide captions for the Movies.

**Movie 1.** Active fluid density from a numerical solution of the dynamic Eqs. (4), (5), and (9) corresponding to Fig. 1(b) for $L = 50$, $\rho_0 = 0.45$, $\tau = 5$, and $\zeta_\rho = 4$. Other parameters as in Table I. Total time is $28 \cdot 10^3$.

**Movie 2.** Active fluid density from a numerical solution of the dynamic Eqs. (4), (5), and (9) corresponding to Fig. 1(c) for $L = 50$, $\rho_0 = 0.6$, $\tau = 5$, and $\zeta_\rho = 4$. Other parameters as in Table I. Total time is $28 \cdot 10^3$.

**Movie 3.** Active fluid density from a numerical solution of the dynamic Eqs. (4), (5), and (9) corresponding to Fig. 1(d) for $L = 50$, $\rho_0 = 0.75$, $\tau = 5$, and $\zeta_\rho = 4$. Other parameters as in Table I. Total time is $28 \cdot 10^3$.

**Movie 4.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for two defects of opposite charges corresponding to Fig. 2(b), $\rho_0 = 0.6$, $\tau = 1$, and $\zeta_\rho = 12$. Other parameters as in Table I. Total time is $2 \cdot 10^3$.

**Movie 5.** Active fluid density and velocity field from a numerical solution of the dynamic Eqs. (4), (5), and (9) for two defects of opposite charges corresponding to Fig. 2(b), $\rho_0 = 0.6$, $\tau = 1$, and $\zeta_\rho = 12$. Other parameters as in Table I. Total time is $2 \cdot 10^3$.

**Movie 6.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for $L = 25$, $\rho_0 = 0.4$, $\tau = 5$, and $\zeta_\rho = 4$. Other parameters as in Table I. Total time is $128 \cdot 10^3$.

**Movie 7.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for $L = 25$, $\rho_0 = 0.5$, $\tau = 5$, and $\zeta_\rho = 4$. Other parameters as in Table I. Total time is $128 \cdot 10^3$.

**Movie 8.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for $L = 25$, $\rho_0 = 0.65$, $\tau = 5$, and $\zeta_\rho = 4$. Other parameters as in Table I. Total time is $128 \cdot 10^3$.

**Movie 9.** Active fluid density and polarity from a nu-

merical solution of the dynamic Eqs. (4), (5), and (9) for $L = 25$, $\rho_0 = 0.7$, $\tau = 5$, and $\zeta_\rho = 4$. Other parameters as in Table I. Total time is $128 \cdot 10^3$.

**Movie 10.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for $L = 25$, $\rho_0 = 0.8$, $\tau = 5$, and $\zeta_\rho = 4$. Other parameters as in Table I. Total time is $128 \cdot 10^3$.

**Movie 11.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for $L = 25$, $\rho_0 = 1$, $\tau = 5$, and $\zeta_\rho = 4$. Other parameters as in Table I. Total time is $128 \cdot 10^3$.

**Movie 12.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for $L = 10$, $\rho_0 = 0.7$, $\tau = 1$, and $\zeta_\rho = 10$. Other parameters as in Table I. Total time is $2 \cdot 10^6$.

**Movie 13.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for $L = 10$, $\rho_0 = 1.3$, $\tau = 0.2$, and $\zeta_\rho = 1$. Other parameters as in Table I. Total time is $2 \cdot 10^6$.

**Movie 14.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for $L = 10$, $\rho_0 = 0.5$, $\tau = 4$, $\zeta_\rho = 4$, $\nu = 0$, and $\zeta_\rho = -1$. Other parameters as in Table I. Total time is $40 \cdot 10^3$.

**Movie 15.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for $L = 10$, $\rho_0 = 0.4$, $\tau = 4$, $\zeta_\rho = 4$, $\nu = -0.5$, and $\zeta_\rho = 1$. Other parameters as in Table I. Total time is $40 \cdot 10^3$.

**Movie 16.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for $L = 10$, $\rho_0 = 0.6$, $\tau = 1$, $\zeta_\rho = 8$, $\nu = 0.5$, and $\zeta_\rho = 0.5$. Other parameters as in Table I. Total time is $40 \cdot 10^3$.

**Movie 17.** Active fluid density and polarity from a numerical solution of the dynamic Eqs. (4), (5), and (9) for $L = 10$, $\rho_0 = 0.4$, $\tau = 2$, $\zeta_\rho = 8$, $\nu = 0$, and $\zeta_\rho = 1$. Other parameters as in Table I. Total time is $40 \cdot 10^3$.

[1] J. Prost, F. Jülicher, and J. F. Joanny, Active gel physics, Nature Physics 11, 111 (2015).

[2] A. B. Verkhovsky, T. M. Svitkina, and G. G. Borisy, Self-polarization and directional motility of cytoplasm., Curr. Biol. 9, 11 (1999).

[3] A. C. Martin, M. Kaschube, and E. F. Wieschaus, Pulsed contractions of an actin-myosin network drive apical constriction, Nature 457, 495 (2009).

[4] M. Mayer, M. Depken, J. S. Bois, F. Jülicher, and S. W. Grill, Anisotropies in cortical tension reveal the physical basis of polarizing cortical flows, Nature 467, 617 (2010).

[5] T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Spontaneous motion in hierarchically assembled active matter, Nature 491, 431 (2012).

[6] G. Duclos, S. Garcia, H. G. Yevick, and P. Silberzan, Perfect nematic order in confined monolayers of spindle-shaped cells, Soft Matter 10, 2346 (2014).

[7] F. Nédélec, T. Surrey, A. C. Maggs, and S. Leibler, Self-

organization of microtubules and motors, Nature 389, 305 (1997).

[8] M. Delarue, J.-F. Joanny, F. Jülicher, and J. Prost, Stress distributions and cell flows in a growing cell aggregate, Interface focus 4, 20140033 (2014).

[9] S. Jain, V. M. Cachoux, G. H. Narayana, S. de Beco, J. D'alessandro, V. Cellerin, T. Chen, M. L. Heuzé, P. Marcq, R.-M. Mège, et al., The role of single-cell mechanical behaviour and polarity in driving collective cell migration, Nature Physics 16, 802 (2020).

[10] J.-M. Armengol-Collado, L. N. Carenza, J. Eckert, D. Krommydas, and L. Giomi, Epithelia are multiscale active liquid crystals, Nature Physics 19, 1773–1779 (2023).

[11] M. Dedenon and K. Kruse, Noise-induced transitions from contractile to extensile active stress in isotropic fluids, Physical Review E 111, 015426 (2025).

[12] T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu,
S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans, and B. Ladoux, Topological defects in epithelia govern cell death and extrusion, *Nature* **544**, 212 (2017).

[13] K. Kawaguchi, R. Kageyama, and M. Sano, Topological defects control collective dynamics in neural progenitor cell cultures, *Nature* **545**, 327 (2017).

[14] Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun, and K. Keren, Topological defects in the nematic order of actin fibres as organization centres of Hydra morphogenesis, *Nature Physics* **17**, 251 (2021).

[15] P. Guillamat, C. Blanch-Mercader, G. Pernollet, K. Kruse, and A. Roux, Integer topological defects organize stresses driving tissue morphogenesis, *Nature Materials* **21**, 588 (2022).

[16] P. G. de Gennes and J. Prost, *The Physics of Liquid Crystals*, 2nd ed., Oxford University Press (Oxford University Press, 2002).

[17] A. J. Vromans and L. Giomi, Orientational properties of nematic disclinations, *Soft Matter* **12**, 6490 (2016), 1507.05588.

[18] X. Tang and J. V. Selinger, Orientation of topological defects in 2D nematic liquid crystals, *Soft Matter* **13**, 5481 (2017), 1706.05065.

[19] D. J. G. Pearce and K. Kruse, Properties of twisted topological defects in 2D nematic liquid crystals†, *Soft Matter* **17**, 7408 (2021), 2104.11293.

[20] C. Blanch-Mercader, P. Guillamat, A. Roux, and K. Kruse, Quantifying Material Properties of Cell Monolayers by Analyzing Integer Topological Defects, *Physical Review Letters* **126**, 028101 (2021).

[21] Y. Maroudas-Sacks, L. Garion, S. Sughanthan, M. Popović, and K. Keren, Confinement Modulates Axial Patterning in Regenerating Hydra, *PRX Life* **2**, 043007 (2024).

[22] Y. Ravichandran, M. Vogg, K. Kruse, D. J. G. Pearce, and A. Roux, Topology changes of Hydra define actin orientation defects as organizers of morphogenesis, *Science Advances* **11**, eadr9855 (2025).

[23] S. P. Thampi, R. Golestanian, and J. M. Yeomans, Velocity Correlations in an Active Nematic, *Phys. Rev. Lett.* **111**, 118101 (2013).

[24] L. Giomi, Geometry and Topology of Turbulence in Active Nematics, *Physical Review X* **5**, 031003 (2015).

[25] R. Alert, J. Casademunt, and J.-F. Joanny, Active Turbulence, *Annual Review of Condensed Matter Physics* **13**, 143 (2022), 2104.02122.

[26] D. J. G. Pearce, Activity Driven Orientational Order in Active Nematic Liquid Crystals on an Anisotropic Substrate, *Physical Review Letters* **122**, 227801 (2019), 1810.04549.

[27] D. J. G. Pearce, Defect order in active nematics on a curved surface, *New Journal of Physics* **22**, 063051 (2020), 2002.06364.

[28] K. Thijssen and A. Doostmohammadi, Binding self-propelled topological defects in active turbulence, *Physical Review Research* **2**, 042008 (2020), 2007.13443.

[29] A. U. Oza and J. Dunkel, Antipolar ordering of topological defects in active liquid crystals, *New Journal of Physics* **18**, 093006 (2016), 1507.01055.

[30] A. Doostmohammadi, M. F. Adamer, S. P. Thampi, and J. M. Yeomans, Stabilization of active matter by flow-vortex lattices and defect ordering, *Nature Communications* **7**, 10557 (2016).

[31] K. Thijssen, M. R. Nejad, and J. M. Yeomans, Role of Friction in Multidefect Ordering, *Physical Review Letters* **125**, 218004 (2020), 2005.01164.

[32] S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti, and V. Vitelli, Topological active matter, *Nature Reviews Physics* **4**, 380 (2022).

[33] T. B. Saw, W. Xi, B. Ladoux, and C. T. Lim, Biological tissues as active nematic liquid crystals, *Advanced materials* **30**, 1802579 (2018).

[34] Y. Maroudas-Sacks and K. Keren, Mechanical patterning in animal morphogenesis, *Annual Review of Cell and Developmental Biology* **37**, 469 (2021).

[35] J. F. Joanny, F. Jülicher, K. Kruse, and J. Prost, Hydrodynamic theory for multi-component active polar gels, *New Journal of Physics* **9**, 422 (2007).

[36] A. C. Callan-Jones and F. Jülicher, Hydrodynamics of active permeating gels, *New Journal of Physics* **13**, 093027 (2011).

[37] J. F. Joanny, K. Kruse, J. Prost, and S. Ramaswamy, The actin cortex as an active wetting layer, *The European Physical Journal E* **36**, 10.1140/epje/i2013-13052-9 (2013).

[38] R. M. Adar and J.-F. Joanny, Active-gel theory for multicellular migration of polar cells in the extra-cellular matrix, *New Journal of Physics* **24**, 073001 (2022).

[39] P. de Gennes and J. Prost, *The Physics of Liquid Crystals*, International Series of Monographs on Physics (Clarendon Press, 1993).

[40] K. Kruse, J.-F. Joanny, F. Jülicher, J. Prost, and K. Sekimoto, Asters, vortices, and rotating spirals in active gels of polar filaments, *Phys. Rev. Lett.* **92**, 078101 (2004).

[41] B. I. Shraiman, Mechanical feedback as a possible regulator of tissue growth, *Proceedings of the National Academy of Sciences* **102**, 3318 (2005).

[42] M. Basan, T. Risler, J.-F. Joanny, X. Sastre-Garau, and J. Prost, Homeostatic competition drives tumor growth and metastasis nucleation, *HFSP journal* **3**, 265 (2009).

[43] S. R. d. Groot and P. Mazur, *Non-Equilibrium Thermodynamics*, Dover Publications Inc (Dover Publications Inc, 1985).

[44] K. Kruse, J.-F. Joanny, F. Jülicher, J. Prost, and K. Sekimoto, Generic theory of active polar gels: a paradigm for cytoskeletal dynamics, *The European Physical Journal E* **16**, 5 (2005).

[45] N. Levernier and K. Kruse, Spontaneous formation of chaotic protrusions in a polymerizing active gel layer, *New Journal of Physics* **22**, 013003 (2020).

[46] L. Giomi, M. J. Bowick, X. Ma, and M. C. Marchetti, Defect Annihilation and Proliferation in Active Nematics, *Phys. Rev. Lett.* **110**, 228101 (2013).

[47] L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek, and M. C. Marchetti, Defect dynamics in active nematics, *Philosophical Transactions of the Royal Society a- Mathematical Physical and Engineering Sciences* **372**, 10.1098/rsta.2013.0365 (2014).

[48] S. Shankar, S. Ramaswamy, M. C. Marchetti, and M. J. Bowick, Defect unbinding in active nematics, *Physical review letters* **121**, 108002 (2018).

[49] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, *Rev Mod Phys* **65**, 851 (1993).

[50] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, *Julia: A fresh approach to numerical computing*, *SIAM review* **59**, 65 (2017).

[51] Code will be made available upon publication.
[52] C. Brito, I. S. Aranson, and H. Chaté, Vortex glass and vortex liquid in oscillatory media, Phys. Rev. Lett. 90, 068301 (2003).

[53] R. Voituriez, J. F. Joanny, and J. Prost, Spontaneous flow transition in active polar gels, EPL (Europhysics Letters) 70, 10.1209/epl/i2004-10501-2 (2007), q-bio/0503022.

[54] G. Duclos, C. Blanch-Mercader, V. Yashunsky, G. Salbreux, J. F. Joanny, J. Prost, and P. Silberzan, Spontaneous shear flow in confined cellular nematics, Nature Physics 14, 728 (2018).

[55] M. Dedenon, C. A. Dessalles, P. Guillamat, A. Roux, K. Kruse, and C. Blanch-Mercader, Density-Polarity Coupling in Confined Active Polar Films: Asters, Spirals, and Biphasic Orientational Phases, Physical Review Letters 131, 268301 (2023).

[56] S. Chandragiri, A. Doostmohammadi, J. M. Yeomans, and S. P. Thampi, Flow States and Transitions of an Active Nematic in a Three-Dimensional Channel, Physical Review Letters 125, 148002 (2020).

[57] I. Lavi, R. Alert, J.-F. Joanny, and J. Casademunt, Dynamical arrest in active nematic turbulence, arXiv 10.48550/arxiv.2407.15149 (2024), 2407.15149.

[58] M. F. Staddon, E. M. Munro, and S. Banerjee, Pulsatile contractions and pattern formation in excitable actomyosin cortex, PLoS Computational Biology 18, e1009981 (2022).

[59] M. Besse, H. Chaté, and A. Solon, Metastability of Constant-Density Flocks, Physical Review Letters 129, 268003 (2022).

[60] B. Lemma, N. P. Mitchell, R. Subramanian, D. J. Needleman, and Z. Dogic, Active Microphase Separation in Mixtures of Microtubules and Tip-Accumulating Molecular Motors, Physical Review X 12, 031006 (2022).

[61] M. Rinaldin, A. Kickuth, B. Dalton, Y. Xu, S. D. Talia, and J. Brugués, Robust cytoplasmic partitioning by solving an intrinsic cytoskeletal instability, bioRxiv, 2024.03.12.584684 (2024).

[62] M. James, D. A. Suchla, J. Dunkel, and M. Wilczek, Emergence and melting of active vortex crystals, Nature Communications 12, 5630 (2021).

[63] D. J. G. Pearce, J. Nambisan, P. W. Ellis, A. Fernandez-Nieves, and L. Giomi, Orientational Correlations in Active and Passive Nematic Defects, Physical Review Letters 127, 197801 (2021).

[64] R. R. Keogh, S. Chandragiri, B. Loewe, T. Ala-Nissila, S. P. Thampi, and T. N. Shendruk, Helical flow states in active nematics, Physical Review E 106, L012602 (2022).

[65] F. Caballero, Z. You, and M. C. Marchetti, Vorticity phase separation and defect lattices in the isotropic phase of active liquid crystals, Soft Matter 19, 7828 (2023).