# Deforming polar active matter in a scalar field gradient

Muhamet Ibrahimi, Matthias Merkel

► To cite this version:

Muhamet Ibrahimi, Matthias Merkel. Deforming polar active matter in a scalar field gradient. New Journal of Physics, 2022, 25 (1), pp.4. 10.1088/1367-2630/acb2e5. hal-03698009v2

HAL Id: hal-03698009

https://hal.science/hal-03698009v2

Submitted on 20 Oct 2023

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>

PAPER • OPEN ACCESS

# Deforming polar active matter in a scalar field gradient

To cite this article: Muhamet Ibrahimi and Matthias Merkel 2023 *New J. Phys.* **25** 013022

View the [article online](#) for updates and enhancements.

You may also like

- [Comparison of contractile and extensile pneumatic artificial muscles](#)

Thomas E Pillsbury, Norman M Wereley and Qinghua Guan

- [A design framework for actively crosslinked filament networks](#)

Sebastian Fürthauer, Daniel J Needleman and Michael J Shelley

- [Defect order in active nematics on a curved surface](#)

D J G Pearce
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED
18 July 2022

REVISED
21 December 2022

ACCEPTED FOR PUBLICATION
13 January 2023

PUBLISHED
26 January 2023

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

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

PAPER

# Deforming polar active matter in a scalar field gradient

Muhamet Ibrahimi ID and Matthias Merkel* ID

Aix Marseille Univ, Université de Toulon, CNRS, CPT (UMR 7332), Turing Centre for Living Systems, Marseille, France

* Author to whom any correspondence should be addressed.

E-mail: matthias.merkel@cnrs.fr

**Keywords:** active matter, morphogenesis, instability, tissue deformation

## Abstract

Active matter with local polar or nematic order is subject to the well-known Simha-Ramaswamy instability. It is so far unclear how, despite this instability, biological tissues can undergo robust active anisotropic deformation during animal morphogenesis. Here we ask under which conditions protein concentration gradients (e.g. morphogen gradients), which are known to control large-scale coordination among cells, can stabilize such deformations. To this end, we study a hydrodynamic model of an active polar material. To account for the effect of the protein gradient, the polar field is coupled to the boundary-provided gradient of a scalar field that also advects with material flows. Focusing on the large system size limit, we show in particular: (a) the system can be stable for an effectively extensile coupling between scalar field gradient and active stresses, i.e. gradient-extensile coupling, while it is always unstable for a gradient-contractile coupling. Intriguingly, there are many systems in the biological literature that are gradient-extensile, while we could not find any that are clearly gradient-contractile. (b) Stability is strongly affected by the way polarity magnitude is controlled. Taken together, our findings, if experimentally confirmed, suggest new developmental principles that are directly rooted in active matter physics.

## 1. Introduction

Active matter is driven out of equilibrium by local injection of mechanical energy, which leads to new properties as compared to inert matter. For instance, active matter with local polar or nematic order is known to exhibit the well-known Simha-Ramaswamy instability [1]. This instability can lead to a spontaneous onset of flows [2] or an instability of the homogeneously deforming state [1, 3], and it has already been observed in several biological systems, including cytoskeletal gels [4], bacterial swarms [5, 6], and cell monolayers *in vitro* [7]. However, it is so far unclear whether this instability appears also *in vivo* during animal morphogenesis, and if not, how it is avoided.

One key process during animal morphogenesis is anisotropic tissue deformation, i.e. pure shear deformation of developing tissue [8–10]. While such deformation can be driven from outside, it is often also driven by active anisotropic stresses generated within the tissue itself [11–18]. To obtain reproducible active anisotropic deformation, rotational symmetry needs to be broken; i.e. there needs to be some kind of directional information encoded in the system, and biological tissues have several ways to do this. For instance, cells in a tissue can possess a polarity, which is defined by an anisotropic distribution of certain polarity proteins within the cell. In developing tissues, such cell polarity often exhibits large-scale ordered patterns [19–21]. In some systems cell polarity can also induce an anisotropic distribution of the motor protein myosin within the cells, and thus control active anisotropic stresses [11, 12]. Such active polar or nematic materials should be prone to the Simha-Ramaswamy instability, and it is so far unclear what prevents it during development.

Animal morphogenesis relies crucially on large-scale protein concentration patterns [8]. Such proteins, called morphogens, are important for the long-range coordination among cells during morphogenesis. In particular, in several tissues, morphogen gradients are also known to control the direction of cell polarity [12, 19, 21–24]. However, it is so far unknown if such protein gradients could help stabilize anisotropic tissue deformation.
From a physics perspective, large-scale dynamics of biological tissues can be described using hydrodynamic active matter models [25, 26]. Such models correspond to expansions describing deviations away from a thermodynamic equilibrium state. They describe materials by averaging over their small-scale features, focusing on the dynamics on long length and time scales. Indeed, in the recent past such hydrodynamic models have been successful in describing the multiple interactions between actively driven tissue deformation, cell polarity, and protein concentration fields [7, 17, 18, 20, 21, 27–34].

Here we use this approach to study under which conditions a gradient of a protein that advects with tissue flows can help stabilize polarity-controlled anisotropic tissue deformation. We describe tissue flows by a velocity field $v$, cell polarity by a polarity field $p$, and the protein concentration field by a scalar field $c$, where we impose a gradient via the boundary conditions. Model details are described in section 2. To examine this model, we first discuss the two limiting cases without polar field (section 3) and without scalar field (section 4), before examining in detail the full model including both scalar and polar fields (section 5). Finally, we discuss under which conditions diffusion and local polarity alignment can stabilize systems of a finite size (section 6).

For the infinite-system-size limit, we find that gradient-contractile systems, i.e. systems with an effective contractile coupling between scalar field gradient and active stresses, are always unstable. Stable tissue deformation is only possible in gradient-extensile systems. Intriguingly, up to one potential exception, effective gradient-extensile coupling is the only coupling that we could find in the biology literature across a multitude of multi-cellular animals. We further show that the stability of deformation strongly depends on how the magnitude of cell polarity is controlled. Taken together, our work suggests new potential developmental principles that are directly rooted in active matter physics.

## 2. Model

### 2.1. Bulk dynamics

We study in 2D the interaction of a scalar field $c(r)$ describing a protein concentration field, a polar field $p(r)$ describing cell polarity, and a velocity field $v(r)$ describing tissue flows (figure 1(a)).

For the scalar field $c$ we focus on simple advective, diffusive dynamics:

$$ \frac{dc}{dt} = D\partial_i^2 c, \qquad (1) $$

where we write the advective time derivative as $dc/dt = \partial c/\partial t + v_i(\partial_i c)$, and $\partial_i$ denotes the partial derivative with respect to the spatial coordinate $r_i$. Here and in the following, we use Einstein notation, and we label spatial dimensions by Latin indices $i, j$.

To describe polarity dynamics, we introduce an effective free energy:

$$ \mathcal{F} = \int \left[ F(p) + \frac{K}{2} (\partial_i p_i)(\partial_i p_i) \right] d^2r, \qquad (2) $$

where the free energy density $F(p)$ controls the polarity magnitude $p = |p|$, and the second term in the integrand controls local polarity alignment, where $K > 0$. The latter is the Frank free energy in the one-constant approximation [35]. Unless stated otherwise, we will set:

$$ F(p) = \frac{A}{2} p^2 + \frac{B}{4} p^4, \qquad (3) $$

where $B > 0$ and $A$ can be either positive or negative. The case $A < 0$ corresponds to a polarity with a finite preferred magnitude (FPM), whereas $A > 0$ corresponds to a polarity with a zero preferred magnitude (ZPM; figure 1(b)).

We use the following polarity dynamics [26, 36]:

$$ \frac{Dp_i}{Dt} = \frac{1}{\gamma}h_i - \nu\tilde{v}_{ij}p_j + \beta\partial_i c. \qquad (4) $$

Here, the time derivative is a co-rotational derivative defined as $Dp_i/Dt = \partial p_i/\partial t + v_j(\partial_j p_i) + \Omega_{ij}p_j$ with the flow vorticity $\Omega_{ij} = (\partial_i v_j - \partial_j v_i)/2$. The first right-hand-side (rhs) term represents the relaxation of the polarity free energy with $\gamma > 0$ being a rotational viscosity and $h_i = -\frac{\delta F}{\delta p_i}$ being the molecular field. The second rhs term is a shear alignment term with coefficient $\nu$ and pure shear rate $\tilde{v}_{ij} = (\partial_i v_j + \partial_j v_i)/2$, where here and in the following, a tilde indicates the symmetric, traceless part of a tensor. For $\nu < 0$, polarity tends to locally align with the extending direction of shear flow, while for $\nu > 0$, polarity tends to locally align with the contracting direction of shear flow. The last rhs term is a coupling to the gradient of the scalar field with
<figure><img src="image_4.png" /></figure>

coefficient $\beta$. Because the dynamics is invariant with respect to the transformation $(\mathbf{p}, \beta) \mapsto (-\mathbf{p}, -\beta)$, we set in the following without loss of generality $\beta > 0$.

We study tissue flow that is governed by incompressible viscous dynamics. We define the stress tensor as

$$ \sigma_{ij} = 2\eta\tilde{v}_{ij} - \Pi\delta_{ij} + \sigma_{ij}^{p} + \tilde{\sigma}_{ij}^{a}, \qquad (5) $$

where $\eta$ is the shear viscosity, $\Pi$ is the hydrostatic pressure, the passive stress due to the polarity $\mathbf{p}$ is

$$ \sigma_{ij}^{p} = \frac{\nu}{2}(p_{i}h_{j} + p_{j}h_{i}) - \frac{1}{2}(p_{i}h_{j} - p_{j}h_{i}), \qquad (6) $$

and we use the following expression for the active stress:

$$ \tilde{\sigma}_{ij}^{a} = \alpha \left( p_{i}p_{j} - \frac{p_{i}^{2}}{2}\delta_{ij} \right). \qquad (7) $$

We use this stress tensor together with force balance and the incompressibility condition:

$$ \partial_{i}\sigma_{ij} = 0 \qquad (8) $$

$$ \partial_{i}v_{i} = 0. \qquad (9) $$

In most of the following, we focus on the infinite-system-size limit, neglecting diffusion $D$ and Frank coefficient $K$. We discuss the effect of finite $D$ and $K$ only in section 6. Moreover, we neglect the passive stress $\sigma_{ij}^{p}$. This corresponds to the limit of a small polarity free energy, $\mathcal{F} \to 0$, while keeping $\mathcal{F}/\gamma$ constant. Otherwise, the passive stress $\sigma_{ij}^{p}$ can be absorbed into a redefinition of $\alpha$ and the hydrostatic pressure.

While more coupling terms between $c$, $\mathbf{p}$, and $\mathbf{v}$ could be added to this model [25, 26], we focus here on the terms directly supported by experimental data on animal morphogenesis [11, 12, 17, 20–22, 32]. Note that the role of the scalar field in this model is different from earlier work, where it represented a concentration of active agents or chemical fuel, and thus scaled the active stresses [25, 37]. Instead, here the scalar field acts as an aligning field for the polarity [1, 3], while also being advected by tissue flows.

## 2.2. Boundary conditions

We study the dynamics in a rectangular periodic box of prescribed time-dependent dimensions $L_x(t) \times L_y(t)$. Because of incompressibility, $L_x(t)L_y(t) = \text{const.}$, and prescribing $L_x(t)$ corresponds to prescribing the average shear rate tensor $\tilde{v}_{ij}^{0}(t)$ with:

$$ \tilde{v}_{xx}^{0}(t) = \frac{1}{L_x} \frac{dL_x}{dt}, \qquad (10) $$

the other diagonal element is $\tilde{v}_{yy}^{0} = -\tilde{v}_{xx}^{0}$, and the off-diagonal elements vanish, $\tilde{v}_{xy}^{0} = \tilde{v}_{yx}^{0} = 0$.
The boundary conditions are periodic for all fields, except for a modification for the scalar field c at the vertical boundary. We set for all $y \in [0, L_y)$:

$$c(0, y) = c(L_x, y) - c_b \qquad (11)$$

with fixed $c_b$. We introduce this modification to ensure that a linear profile $c = xc_b/L_x$ is stationary.

When prescribing not the box dimensions but an external stress anisotropy $\tilde{\sigma}_{ij}^{\text{ext}}(t)$, equation (5) implies that the system will shear with rate $\tilde{v}_{xx}^0(t) = (\tilde{\sigma}_{xx}^{\text{ext}}(t) - \langle \tilde{\sigma}_{xx}^a \rangle)/2\eta$. We use $\tilde{\sigma}_{ij}^{\text{ext}} = \langle \tilde{\sigma}_{ij} \rangle$, where $\tilde{\sigma}_{ij}$ is the symmetric, traceless part of $\sigma_{ij}$, and $\langle \cdot \rangle$ is the spatial average over the system. In the following, we will study stability for any constant $\tilde{v}_{xx}^0(t)$. Two specific cases are: (a) a system with fixed size, $\tilde{v}_{xx}^0 = 0$, and (b) a freely deforming system, $\tilde{\sigma}_{ij}^{\text{ext}}(t) = 0$, which has the box deformation rate $\tilde{v}_{xx}^0(t) = -\langle \tilde{\sigma}_{xx}^a \rangle/2\eta$.

## 2.3. Dimensionless dynamics

We non-dimensionalize the dynamics, equations (1)–(11), by choosing $c_b$ as unit for the scalar field, $\sqrt{|A|/B}$ as polarity unit, $L_x(0)$ as length scale, $|\alpha A|/B$ as stress scale, and $\eta B/|\alpha A|$ as time scale. In the rest of this article, we use the accordingly rescaled dimensionless quantities. The dimensionless dynamical equations are

$$\frac{\mathrm{d}c}{\mathrm{d}t} = D\partial_i^2 c \qquad (12)$$

$$\frac{\mathrm{D}p_i}{\mathrm{D}t} = -\frac{g(p)}{\tau}p_i - \nu\tilde{v}_{ij}p_j + \beta\partial_i c + \kappa\partial_j^2 p_i \qquad (13)$$

$$0 = \partial_i^2 v_j - \partial_j \Pi' + \operatorname{sgn}(\alpha) \partial_i (p_i p_j). \qquad (14)$$

$$\partial_i v_i = 0. \qquad (15)$$

Here, we defined $g(p) = F'(p)/|A|p$ such that with our choice in equation (3):

$$g(p) = \operatorname{sgn}(A) + p^2. \qquad (16)$$

Here we have introduced the sign function, $\operatorname{sgn}(A) := A/|A|$. We moreover introduced the dimensionless time scale $\tau = \gamma/|A|$ over which the polarity magnitude relaxes, a polarity alignment coefficient $\kappa = K/\gamma$, and we set $\Pi' = \Pi + p^2/2$.

In dimensionless units, the modified boundary condition becomes $c(0, y) = c(L_x, y) - 1$, the box deformation rate for the freely deforming system is $\tilde{v}_{xx}^0 = -\operatorname{sgn}(\alpha)\langle p_x^2 - p_y^2 \rangle/4$, and we define the box shear $l_x(t) := L_x(t)/L_x(0) \equiv L_x(t)$.

## 3. Scalar field only

We first discuss the special case without polar field. This can also be regarded as the limit where the polarity relaxes to the scalar field gradient adiabatically fast:

$$p_i = \partial_i c. \qquad (17)$$

More precisely, this corresponds to the limit $\beta \to \infty$ while keeping $\beta\tau = 1$, and using the polarity potential $F(p) = p^2/2$ (i.e. $g(p) = 1$). As a consequence, the active stress is given by the gradient of the scalar field $\tilde{\sigma}_{ij}^a = \operatorname{sgn}(\alpha) [(\partial_i c)(\partial_j c) - (\partial_k c)^2 \delta_{ij}/2]$, like in Active Model H [38].

Example numerical solutions of the dynamics for the freely deforming system are shown in figure 2(a): gradient-contractile systems, $\alpha > 0$, are unstable, while gradient-extensile systems, $\alpha < 0$, are marginally stable.

### 3.1. Fixed system size, $\tilde{v}_{xx}^0 = 0$

We first perform a linear stability analysis for a non-deforming system, $\tilde{v}_{xx}^0 = 0$, by perturbing the scalar field c around the linear profile:

$$c = c_0 + \delta c, \qquad (18)$$

where we defined the linear profile $c_0$ as:

$$c_0(\mathbf{r}, t) = \frac{x}{l_x(t)} \qquad (19)$$
<figure><img src="image_5.png" /></figure>

with $\mathbf{r} = (x, y)$. We obtain for the growth rate $\omega$ of a perturbation with wave vector $\mathbf{k} = k(\cos \phi, \sin \phi)$ (figure 2(d), appendix B.2):

$$ \omega(\mathbf{k}) := \frac{\partial_t \delta c(\mathbf{k}, t)}{\delta c(\mathbf{k}, t)} = \operatorname{sgn}(\alpha) \sin^2 \phi. \qquad (20) $$

Thus, consistent with earlier work [39], the system is linearly unstable for gradient-contractile systems ($\alpha > 0$), while it is marginally stable for gradient-extensile systems ($\alpha < 0$), where the only modes that do not decay over time are those with wave vector angles $\phi = 0$ and $\phi = \pi$. The magnitudes of growth and decay rates can thus be up to $\pm 1$ in dimensionless units, which is four times the magnitude of the free deformation rate, which is $\tilde{v}_{xx}^{\text{free}} = -\operatorname{sgn}(\alpha)/4$.

Intuitively, the difference in behavior between gradient-extensile and gradient-contractile cases is illustrated in figures 2(b) and (c), which represent the behavior of a perturbation with wave vector angles $\phi = \pi/2$. A perturbation $\delta c$ of the scalar field (red symbols in figures 2(b) and (c)) induces to first order a perturbation in the active stress nematic of $\delta \tilde{\sigma}_{xi}^a \sim k_i \delta c$, which have angles of $\pm \pi/4$ (black double arrows). For the gradient-extensile case ($\alpha < 0$), figure 2(c)), this stress perturbation generates flows (blue arrows) that advect regions with positive $\delta c$ in positive $x$ direction, leading to a local decrease in $c$ due to the overall $c$ gradient. Hence, for the gradient-extensile case this perturbation decays. Analogously, for the gradient-contractile case, advection leads to an amplification of the $\pi/2$ mode (figure 2(b)).

### 3.2. Deforming system, $\tilde{v}_{xx}^0 \neq 0$

For a deforming system, the linear profile $c_0$ (equation (19)) becomes distorted by the system's overall deformation, which means that it is not a stationary state with respect to the lab frame anymore. However, stationarity is still possible with respect to 'co-deforming' coordinates $\tilde{\mathbf{r}} = (\tilde{x}, \tilde{y})$, which we define as (see appendix A):
$$ \bar{x} = \frac{xL_x(0)}{L_x(t)} \equiv l_x^{-1}(t)x \qquad (21) $$

$$ \bar{y} = \frac{yL_y(0)}{L_y(t)} \equiv l_y(t)y. \qquad (22) $$

These coordinates thus map any point $r$ in the system at time $t$ to its affinely rescaled position $\bar{r}$ at time zero. As a consequence, the linear profile $c_0$ does not become distorted when written in co-deforming coordinates

$$ c_0(\bar{r}, t) = \bar{x}. \qquad (23) $$

Hence, in co-deforming coordinates $c_0$ can be regarded as stationary.

Co-deforming coordinates are useful also because they help us solve the linearized dynamics of the system, which is required to analyze the system's stability (appendix A.3). We find that the solutions of the linearized dynamics are 'co-deforming Fourier modes':

$$ \delta c(\bar{r}, t) = \int \delta c(\bar{k}, t) e^{i\bar{k}\bar{r}} d^2\bar{k}. \qquad (24) $$

We call $\bar{k}$ a co-deforming wave vector. A co-deforming Fourier mode with wave vector $\bar{k}$ is distorted over time by the overall system deformation (figure 3(d), appendix A.1):

$$ k_x(\bar{k}, t) = l_x^{-1}(t)\bar{k}_x \qquad (25) $$

$$ k_y(\bar{k}, t) = l_y(t)\bar{k}_y. \qquad (26) $$

These relations will introduce a time dependence in the growth rate of co-deforming Fourier modes.

For the deforming scalar-only system, we obtain for the growth rate of a co-deforming Fourier mode with wave vector $\bar{k}$ (see appendix B.2):

$$ \omega(\bar{k}, t) := \frac{\partial_t \delta c(\bar{k}, t)}{\delta c(\bar{k}, t)} = \operatorname{sgn}(\alpha) \frac{\sin^2 \phi(\bar{k}, t)}{l_x^2(t)}. \qquad (27) $$

Here, the angle $\phi$ of the lab-frame wave vector $\bar{k}$ depends on $\bar{k}$ and $t$ through equations (25) and (26). As a consequence of the time-dependent right-hand side of equation (27), the solutions $\delta c(\bar{k}, t)$ are generally not exponential in $t$ any more. However, we can still discuss the stability of the system. Indeed, like in the case with fixed system size, the system is unstable in the gradient-contractile case, $\alpha > 0$, while it is marginally stable in the gradient-extensile case, $\alpha < 0$.

## 4. Polar field only

We now revisit the limit where the scalar field plays no role, $\beta = 0$, i.e. where our systems becomes a polar active matter system. While such systems have been discussed before in the literature [1, 2, 25], we will shed light on a few new aspects.

Example simulations of a freely deforming system with fixed polarity magnitude are shown in figure 3(a): we find unstable behavior independent of the sign of the shear alignment $\nu$.

In our discussion below, we focus here on the extensile case, $\alpha < 0$. For $\beta = 0$ these results can be directly mapped to the contractile case.

### 4.1. Fixed system size, $\tilde{v}_{xx}^0 = 0$

Here we focus exclusively on the case of a finite preferred polarity magnitude (FPM, figure 1(b)). For fixed system size, $\tilde{v}_{xx}^0 = 0$, there are no stationary solutions with finite polar order for ZPM.

We briefly revisit the stability of the polar ordered state [1, 2]:

$$ \mathbf{p}_0 = \dot{\mathbf{x}}, \qquad (28) $$

$$ \mathbf{v}_0 = 0. \qquad (29) $$

The system is unstable for every value of the shear alignment coefficient $\nu$ [1]. To prepare for the later sections, we now revisit the intuitive explanation of this for the limit of $\tau \to 0$, where the polarity magnitude is fixed to one, $|\mathbf{p}| = 1$ (figures 3(b) and (f), appendix B.3.1).
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Linear stability of the polar-only system, $\beta = 0$. (a) Simulations of a freely deforming system with fixed-magnitude polarity and extensile active stresses, for two different values of the shear alignment parameter $\nu$. (b) Growth rate $\omega$ of a perturbation with angle $\phi$ in a system with a fixed polarity magnitude (i.e. sgn($A$) = -1 and $\tau \to 0$) for three different $\nu$ values. (c)–(e) Schematics illustrating the perturbation dynamics for three different $\phi$ (section 4.1). (f) Perturbation growth rate maximized across all wave vectors, $\omega_{max}$, for fixed system (blue) and freely deforming system (purple). The dashed lines represent the growth rates of the respective bend modes ($\phi = 0$). Green stars correspond to the $\nu$ values used for the simulations in panel a. (g) Illustration of the change of the wave vector $\mathbf{k}$ measured in the lab frame (arrow) for a given co-deforming mode $\tilde{\mathbf{k}}$ as the system is sheared from $l_x(t) = 1$ (A) to $l_x(t) = 2$ (B; compare equations (25) and (26)). (h) Perturbation angle $\phi$ measured in the lab frame versus the box shear $l_x$, for three different co-deforming angles $\tilde{\phi}$. (i) Change of the growth rates $\omega$ of three co-deforming perturbation modes (same as in panel h) as the box width $l_x$ changes during the box deformation. The co-deforming $\tilde{\phi} = \pi/4$ mode transiently experiences a positive growth rate (positive $\omega$ values between '+' signs). Fixed polarity magnitude, $\tilde{v}_{xx}^0 = 1/8$, $\nu = 3$. (j)–(m) Stability phase diagrams for (j), (k) fixed polarity magnitude and (l), (m) zero preferred polarity magnitude with $\tau = 10$, where red, gray, and yellow respectively indicate unstable regions, stable regions, and stable regions with transiently growing modes. The green dashed line in (j) represents the deformation rate of the freely deforming system. Black lines marking phase boundaries are analytical curves derived in appendix B.3. (k)–(l) Maximal perturbation amplification factor $S$ in the stable regions with transiently growing modes. Orange star in (k) corresponds to the parameters used for panel i.</figcaption></figure>
For $\nu < 1$ the bend mode $\phi = 0$ is unstable (figure 3(b) and blue dashed line in figure 3(f)) [2]. The mechanism driving this instability is illustrated in figure 3(c): For fixed polarity norm, perturbations in the polarity are exclusively orientational, i.e. $\delta p$ is oriented along the y axis (green arrows). This affects the linear perturbation in the active stress tensor $\delta\tilde{\sigma}_{ij}^a$, which is in general given by:

$$ \delta\tilde{\sigma}_{xi}^a = \alpha p_0 \delta p_i. \qquad (30) $$

For $\delta p$ oriented along the y axis, equation (30) implies that $\delta\tilde{\sigma}_{ij}^a$ is oriented at angles of $\pm\pi/4$ with respect to the x axis (black arrows). This active stress perturbation creates flow (blue arrows), which for $\nu < 1$ acts through a combination of co-rotational and shear alignment effects to amplify the perturbation $\delta p_y$.

For $\nu \ge 1$, where the bend mode is stable, there are always other unstable modes [1]. To see this, we first discuss the mode with $\phi = \pi/4$ (figure 3(d)). In this case, the angle of the wave vector $k$ is parallel or perpendicular to the angle of the active stress perturbation $\delta\tilde{\sigma}_{ij}^a$, and thus no flow appears due to incompressibility. However, when slightly decreasing $\phi$ below $\pi/4$ (figure 3(e)), flow appears, which acts through the co-rotational effect to amplify $\delta p_y$. This co-rotational effect always dominates over the shear alignment when approaching $\phi = \pi/4$. As a consequence, there are always perturbations with angles $\phi$ below $\pi/4$ whose growth rate are positive.

These results do not fundamentally change when allowing for a finite value of the polarity magnitude relaxation time, $\tau > 0$, i.e. a ‘soft’ polarity magnitude. For $\nu < 1$ the fastest growing mode is still the bend mode, whose growth rate is unaffected by $\tau$. For $\nu > 1$, the growth rate of the fastest growing mode decreases with increasing $\tau$. However it is still positive for any $\tau$ (figure B7). Hence, the system with fixed size is always unstable, also for finite $\tau$.

## 4.2. Deforming system, $\tilde{v}_{xx}^0 \neq 0$

For a deforming system with given box shear rate $\tilde{v}_{xx}^0$ the stationary, homogeneously deforming state is:

$$ p_0 = p_0 \hat{x}, \qquad (31) $$

$$ v_{0,i} = \tilde{v}_{ij}^0 r_j, \qquad (32) $$

where the value of the stationary polarity magnitude $p_0$ depends on the given box shear rate $\tilde{v}_{xx}^0$. In this section we discuss the stability of this stationary state using co-deforming perturbations (section 3.2, appendix A).

### 4.2.1. Finite preferred polarity magnitude (FPM)

For simplicity, we focus here on fixed-magnitude polarity ($\tau \to 0$). For fixed box shear rate $\tilde{v}_{xx}^0$, the growth rate of a solution of the linearized dynamics with co-deforming wave vector $\bar{k}$ is (appendix B.3.1):

$$ \frac{\partial_t \delta p_y(\bar{k}, t)}{\delta p_y(\bar{k}, t)} = \omega(\phi) \quad \text{with} \quad \omega(\phi) := -\frac{1}{2} \cos(2\phi) [\nu \cos(2\phi) - 1] + 2\nu \tilde{v}_{xx}^0. \qquad (33) $$

Here, $\phi = \phi(\bar{k}, t)$ is the angle of the wave vector measured in the lab frame, $\bar{k}$, which varies with time $t$ according to equations (25) and (26) (figures 3(g) and (h)). Only bend and splay modes, $\bar{\phi} \in \{0, \pi/2, \pi, 3\pi/2\}$, have angles $\phi$ that are independent of time: $\phi = \bar{\phi}$ (figure 3(h)).

The linear stability phase diagram according to equation (33) is shown in figure 3(j) depending on box shear rate $\tilde{v}_{xx}^0$ and shear alignment $\nu$. Red regions indicate unstable systems, such as systems with fixed size, $\tilde{v}_{xx}^0 = 0$, as we have seen in the previous section, and freely deforming systems, $\tilde{v}_{xx}^0 = 1/4$ (green dashed line). Dark gray regions indicate stable systems, which occur for deformation rates $\tilde{v}_{xx}^0$ larger than the free deformation or for negative deformation rates $\tilde{v}_{xx}^0 < 0$. The latter corresponds to the case where the boundary conditions force the system to deform perpendicular to the free deformation.

In the yellow parameter regions in figure 3(j), the system is stable except for only transiently growing modes. In these regions, both bend and splay modes are decaying, while a co-deforming perturbation mode with given $\bar{k}$ that is neither bend nor splay mode can transiently experience a positive growth rate. This happens as the angle of the corresponding wave vector measured in the lab frame, $\phi(\bar{k}, t)$, which varies as the box is sheared, passes through a regime of angles with a positive growth rate $\omega(\phi)$ (figures 3(b) and (g)-(i)). For positive (negative) $\tilde{v}_{xx}^0$, all angles $\phi(\bar{k}, t)$ ultimately approach splay (bend) modes for $t \to \infty$, which have negative growth rates. Hence, the mode grows only transiently.
The amplitude of a transiently growing co-deforming mode will only increase by some maximal amplification factor until it decreases again (figure 3(i)). This factor attains the same maximal value S for all co-deforming modes $\bar{k}$ that pass both positive zeros in the $\omega(\phi)$ curve (figure 3(b)):

$$S = \exp \left[ \int_{t_1}^{t_2} \omega(\phi(\bar{\phi}, t)) dt \right], \qquad (34)$$

where $t_1$ and $t_2$ are the times when $\omega(\phi(\bar{\phi}, t))$ passes zero (marked by ‘+’ signs in figure 3(i)). In figure 3(k), we plot the maximal amplification factor depending on $\nu$ and $\tilde{v}_{xx}^0$. This amplification factor is affected by two parameters, the area of the positive region in $\omega(\phi)$ (figure 3(b)) and the speed by which it is traversed, which is set by the box shear rate $\tilde{v}_{xx}^0$ (figures 3(g) and (h)). The latter is the reason that we see a diverging amplification factor as $\tilde{v}_{xx}^0 \to 0$ (figure 3(k)).

### 4.2.2. Zero preferred polarity magnitude (ZPM)

For ZPM, there is no stationary state with finite polar order in the regime where $\nu\tilde{v}_{xx}^0\tau \ge -1$ (white region in figure 3(m). Stabilizing polarity in these regions would require higher-order terms in the polarity free energy [40], which we neglect for brevity here. We focus here on the parts of the parameter space where a stationary state finite polarity magnitude $p_0$ exists, where $p_0 = \sqrt{-1 - \nu\tilde{v}_{xx}^0\tau}$.

We find that the linear stability diagram for systems with ZPM is quite different from systems with a fixed polarity norm (compare figures 3(j) and (m)). For instance, large parts of the regime with $\nu < 0$ are unstable for ZPM, which is different from the case with a fixed polarity norm, where the system is stable in this regime. This difference comes mostly from the fact that the polarity norm $p_0$ can become much larger than one for ZPM (appendix B.3.2). As a consequence, the perturbation growth rate is dominated by the flow created by the active stress, which scales as $\sim p_0^2$ and destabilizes polarity in this regime (appendix B.3.2). Conversely, it can be shown that for $\nu \ge 1$ the system is always stable (appendix B.3.2). Taken together, in the regime of zero preferred polarity magnitude, the linear stability phase diagram is significantly affected by the fact that the polarity magnitude $p_0$ is influenced by the box shear rate through shear alignment.

### 4.3. Comparison to system with scalar field only

In polar-only systems (this section, 4), the system's behavior is symmetric with respect to whether the active stress is extensile or contractile. However, this was not true in scalar-only systems (section 3) [38, 39]. Where does this difference come from?

To address this question, we map the scalar field gradient to an effective polar field $q$ with $q_i = \partial_i c$. This transforms the scalar field dynamics, equation (12) into:

$$\frac{Dq_i}{Dt} = -\tilde{v}_{ij}q_j, \qquad (35)$$

and the active stress becomes $\tilde{\sigma}_{ij}^a = \alpha(q_i q_j - q^2 \delta_{ij}/2)$. Thus, for an incompressible system, the scalar field dynamics corresponds to the dynamics of a polar field in the limit of no magnitude control, $\tau \to \infty$, and with a shear alignment coefficient of +1. Indeed, using our fixed-system-size results for the polar system in this limit (figure B7), we can retrieve our stability results for both extensile and contractile cases of the scalar system. Hence, it is the effective positive shear alignment coefficient of +1 that breaks the extensile/contractile symmetry in the scalar-only system.

## 5. Scalar and polar field

Here we examine the general case where the interactions between scalar, polar, and flow fields play a role. We discuss the stability of the homogeneously deforming state given by

$$c_0 = \frac{x}{l_x(t)}, \qquad (36)$$

$$\mathbf{p}_0 = p_0 \hat{\mathbf{x}}, \qquad (37)$$

$$v_{0,i} = \tilde{v}_{ij}^0 r_j. \qquad (38)$$

In the full system, the polarity magnitude $p_0$ generally depends on time, because it is coupled to the scalar field gradient, which is constantly becoming flatter. Thus, to simplify our discussion in this section, we choose a time-dependent $\beta$ that compensates for the flattening gradient:

$$\beta(t) = \beta_0 l_x(t). \qquad (39)$$
<figure><img src="image_7.png" /><figcaption>Figure 4. Linear stability of the non-deforming full system. A coupling to a scalar field gradient can suppress the instability of an active polar system. This stabilization is more efficient for ZPM than for FPM. (a) Result summary: For contractile active stresses, the system is always unstable (left). However, for extensile active stresses, the presence of the scalar field can stabilize the system (right). (b) Stationary polarity magnitude $p_0$ vs. $\beta_0\tau$ shown for finite (dark green) and zero (light green) preferred polarity magnitude. (c) Stiffness ratio G vs. $\beta_0\tau$ shown for both types of polarity. (d), (e) Maximal perturbation growth rate $\omega_{max}$ shown depending on coupling to the scalar field gradient, $\beta_0$, and shear alignment, $\nu$, for (d) fixed polarity magnitude ($\tau \to 0$) and (e) zero preferred polarity magnitude with $\beta_0\tau = 1$. The system with fixed polarity magnitude is always unstable (red region), even for strong coupling to the scalar field gradient, whereas the system with zero preferred polarity magnitude can be marginally stabilized (light gray regions) by a strong enough coupling to the scalar field gradient $\beta_0$. The black solid curves indicate the places where the bend mode growth rates cross zero. (f) Perturbation growth rate $\omega/p_0^2$ over wave vector angle $\phi$ in the limit of strong coupling to the scalar field, $\beta_0 \gg p_0^3$, shown for finite (dark green) and zero (light green) preferred polarity magnitude. Solid curves: $\beta_0\tau = 10^{-2}$, dashed curves: $\beta_0\tau = 5$. (g) Stress perturbation angle $\theta_\sigma$ over wave vector angle $\phi$, shown for the same cases discussed in panel e. (h)–(k) Sketches illustrating the linear stability in the limit of strong coupling to the scalar field, $\beta_0 \gg p_0^3$, for finite (h), (i) and zero (j), (k) preferred polarity magnitude. Perturbation angles are $0 < \phi < \pi/4$ in (h), (j) and $\pi/4 < \phi < \pi/2$ in (i), (k).</figcaption></figure>

In this case, the state defined by equations (36)–(38) is stationary in co-deforming coordinates, where $p_0$ is implicitly set by:

$$g(p_0) = \frac{\beta_0\tau}{p_0} - \nu\tilde{v}_{xx}\tau. \qquad (40)$$

We prove in appendix B.1 that gradient-contractile systems, i.e. systems with $\beta > 0$ and $\alpha > 0$, are always unstable (figure 4(a) left). In the present section, we thus focus our discussion on the gradient-extensile case, $\beta > 0$ and $\alpha < 0$ (figure 4(a) right), where we show that the scalar gradient can indeed stabilize the Simha-Ramaswamy instability under certain conditions. Such a stabilization depends in particular on how polarity magnitude is controlled; it is more effective for ZPM than for FPM.
## 5.1. Fixed system size, $\tilde{v}_{xx}^0 = 0$

### 5.1.1. Finite preferred polarity magnitude (FPM)

We first discuss the limit of a fixed polarity norm, $\tau \to 0$, where $p_0 = 1$. In this case, the system behavior depends only on two dimensionless parameters, $\nu$ and $\beta_0$. We find that in this case the system is always unstable (figure 4(d) left), which we explain in the following.

We discuss two limits, weak and strong coupling of polarity to the scalar field gradient, $\beta_0$. Weak coupling, $\beta_0 \ll p_0^3$ (appendix B.4.1), corresponds to the polarity-only case, which we discussed in section 4. In this case, polarity was always unstable for fixed system dimensions (compare figure 4(d) right to figure 3(d)).

In the limit of strong coupling to the scalar field gradient, $\beta_0 \gg p_0^3$, one might expect the same result as for the scalar-only system, section 3, which was marginally stable in the gradient-extensile case. However, we find that this is not the case, and the full system is unstable instead (figure 4(d) right). To intuitively understand why, consider a given perturbation $\delta c$ of the concentration field along any wave vector with angle $0 < \phi < \pi/4$ (red symbols in figure 4(h)). The polarity perturbation $\delta p$ (green arrows) will adjust to $\delta c$ adiabatically quickly for strong coupling $\beta_0$. However, because of the fixed polarity magnitude, $\delta p$ will point in $\pm \hat{y}$ direction. As a consequence, the perturbation in the active stress tensor $\delta \tilde{\sigma}_{ij}^a$ is oriented along angles of $\pm \pi/4$ (black arrows, using equation (30)). The resulting flow (blue arrows) has a component that points in $-\hat{x}$ direction in regions where $\delta c$ is positive. Thus, due to convection and the gradient in $c_0(x)$, the amplitude of $\delta c$ increases. Hence, the system is unstable for fixed polarity magnitude.

These ideas generalize to the case of a finite polarity relaxation time $\tau$. For given scalar field perturbation $\delta c$, and strong coupling to the scalar field, $\beta_0 \gg p_0^3$, the polarity perturbation $\delta p$ relaxes adiabatically to (appendix B.4.1):

$$ \delta p_x = iGp_0k_x\delta c \quad (41) $$

$$ \delta p_y = ip_0k_y\delta c \quad (42) $$

with

$$ G = \frac{g(p_0)}{g(p_0) + p_0 g'(p_0)}. \quad (43) $$

Thus, the polarity perturbation $\delta p$ does not locally align parallel to the gradient of $\delta c$. Instead, the angle $\theta_p$ of the polarity perturbation is given by (precise definition of $\theta_p$ in appendix B.4.1):

$$ \tan \theta_p = \frac{1}{G} \tan \phi. \quad (44) $$

The prefactor $G$ arises in equations (41) and (42), because for a finite preferred polarity norm, the polarity free energy $F(p)$ panelizes perturbations $\delta p_x$ stronger than perturbations $\delta p_y$. In other words, $G$ is the ratio between the effective stiffnesses associated with changes in polarity away from the stationary state along $y$ and $x$ axes. Using equation (40) with $\tilde{v}_{xx}^0 = 0$, the value of $G$ depends on the product $\beta_0\tau$ only (figure 4(c)); for finite preferred polarity magnitude, $G$ increases from $G=0$ at $\beta_0\tau \to 0$ to maximally $G \to 1/3$ at $\beta_0\tau \to \infty$.

The stiffness ratio $G$ controls the stability of the system. To see this, we discuss the flow created by the polarity perturbation $\delta p$, whose $x$ component is given by (appendix B.4.1):

$$ \delta v_x = -p_0^2 \hat{G} \sin \phi \sin(2[\theta_\sigma - \phi]) \delta c. \quad (45) $$

Here, $\hat{G} > 0$, and $\theta_\sigma$ is the angle of the active stress perturbation nematic $\delta \tilde{\sigma}_{ij}^a$ (precise definitions in appendix B.4.1). Thus, the relative direction between active stress perturbation angle $\theta_\sigma$ and wave vector angle $\phi$ determines the direction of the flow in $x$ direction. Because the flow in equation (45) advects the scalar field, the system is unstable whenever $\delta v_x/\delta c < 0$. Thus, equation (45) implies that the system is unstable whenever $\phi < \theta_\sigma$, where $\theta_\sigma$ can be obtained using equation (30) as: $\theta_\sigma = \theta_p/2$.

For example for fixed polarity magnitude, i.e. $\beta_0\tau \to 0$, we have $G=0$ (figure 4(c)), which implies with equation (44) that $\theta_p = \pi/2$. Further, $\theta_\sigma = \theta_p/2 = \pi/4$, and thus the modes with $\phi < \theta_\sigma = \pi/4$ are unstable (figures 4(f) and (g)).

Also for finite $\beta_0\tau$ there will always be angles $\phi$ for which $\phi < \theta_\sigma$. This is because for $\phi \ll 1$, equation (44) implies that $\theta_\sigma = \theta_p/2 \simeq \phi/(2G)$. Since we have $G < 1/3$ (figure 4(c), appendix B.4.1), it directly follows that $\phi < \phi/(2G) = \theta_\sigma$ for small $\phi$ (see dashed dark curve in figure 4(g)), and thus the system is unstable.

Taken together, for finite preferred polarity magnitude (FPM), the system is always unstable, even in the limit of strong coupling to the scalar field gradient. This is ultimately because $G < 1/3$, i.e. because polarity perturbation is less than 1/3 as stiff along the $y$ axis than along the $x$ axis.
<figure><img src="image_8.png" /></figure>

**Figure 5.** Linear stability of the deforming full system. Stability phase diagrams depending on shear alignment $\nu$ and box shear rate $\tilde{v}_{xx}^0$. (a)–(c) Fixed polarity magnitude ($\tau \to 0$) for $\beta_0 = (0.1, 1, 10)$. As the coupling $\beta_0$ to the scalar field gradient increases, the extent of unstable regions (red) shrink. (d)–(f) Zero preferred polarity magnitude ($\beta_0\tau = 1$). As $\beta_0$ increases, the extent of the marginally stable region (light gray) increases. The dashed green curves indicate the deformation rate of the freely deforming system. In panel d, this line is undefined for some negative $\nu$, because polarity for the freely deforming system would diverge there. Phase boundaries displayed as solid black curves are derived in appendix B.4.2. Note the different scaling of the color map between panels (a)–(c) vs. (d)–(f).

### 5.1.2. Zero preferred polarity magnitude (ZPM)

For ZPM and weak coupling to the scalar field, $\beta_0 \ll p_0^3$, there is no polar-only case that we can directly compare to, because the stationary state has zero polarity in this case. However, we find that the system is unstable for negative $\nu$ (figure 4(e)), which is due to an unstable bend mode (indicated by the solid black curve).

In the limit of strong coupling to the scalar field gradient, $\beta_0 \gg p_0^3$, stability is again controlled by the stiffness ratio $G$. However, because for ZPM the stiffness ratio $G$ is more isotropic, attaining values closer to one (figure 4(c)), the system can become marginally stable (figures 4(j) and (k)). This can be seen following the same line of argument as for FPM. In particular, it can be shown that the system is marginally stable whenever $G \ge 1/2$, i.e. whenever polarity perturbation is at least $1/2$ as stiff along the $y$ axis than along the $x$ axis (see appendix B.4.1). This is the case whenever $\beta_0\tau \le 2$ (figures 4(c) and (f)), i.e. when polarity magnitude relaxation is fast enough as compared to the effect of the scalar field gradient.

## 5.2. Deforming system, $\tilde{v}_{xx}^0 \ne 0$

### 5.2.1. Finite preferred polarity magnitude (FPM)

Here we focus on the case of fixed polarity norm, $\tau = 0$. In this case, the polar-only system had two stable regions for $\nu \tilde{v}_{xx}^0 < 0$ (figure 3(j)). However, when adding a weak coupling to the scalar field, $\beta_0 \ll p_0^3$, these previously stable regions now turn into regions with transiently growing modes (figure 5(a)). This is because when including the scalar field, an additional $\omega(\phi)$ branch appears, and thus more modes that could potentially grow. In particular, while perturbations in the polar field $\delta p_y$ relax to zero relatively quickly, any perturbation in the scalar field $\delta c$ induces a small perturbation in the polarity of $\delta p_y = -ik_y\beta_0\delta c/\omega$. This polarity creates flows that then advect the scalar field. As a consequence, a scenario qualitatively similar to the one discussed in figure 4(h) arises. Here, this gives rise to growing modes with $0 < \phi < \pi/4$ that are only transiently growing, because the system is deforming. Their amplification factor is small, $S \approx 1$, because the magnitude of the growth rate is proportional to $\beta_0$ in this case.
When increasing the coupling $\beta_0$, the region of parameter space where the system is stable with transiently growing modes is expanding (see figures 5(b) and (c)). In the limit $\beta_0 \gg p_0^3$, modes with $0 < \phi < \pi/4$ are again unstable for the same reason as discussed in section 5.1.1 (figure 4(g)). In this limit, we have transiently growing modes, whose amplification factor $S$ is independent of shear alignment $\nu$ and scales as $S \sim 1/|\tilde{v}_{xx}|^2$.

### 5.2.2. Zero preferred polarity magnitude (ZPM)

For ZPM and weak coupling to the scalar field gradient, $\beta_0 \ll p_0^3$, stability is somewhat similar to the polarity-only case (compare figures 3(m) and 5(d)). The main difference is that most of the regions of parameter space without stationary state in the polarity-only system become unstable when adding a weak coupling to a scalar field gradient. Increasing the coupling to the scalar field $\beta_0$ generally leads to an expansion of the parameter regime with marginal stability (see figures 5(e) and (f)).

## 6. System size effects

Until now we have discussed the limit of infinite system size, where neither diffusion $D$ nor polarity alignment $\kappa = K/\gamma$ can suppress large-wavelength perturbations with $k \to 0$. However, a finite system size limits the wavelength from above, which can allow diffusion and polarity alignment to stabilize the system. Here we discuss under which conditions they can do so, focusing on the case of a fixed system size and a fixed polarity magnitude. In this case, the system was always unstable in the infinite-system-size limit (section 5.1.1).

For the limits of scalar-only and polar-only dynamics, the stability criteria are already known. For the scalar-only limit, an otherwise unstable gradient-contractile system is stabilized whenever [39]

$$Dk^2 > 1, \qquad (46)$$

where we use dimensionless units (section 2.3), i.e. rates are in units of $|\alpha|/\eta$. For the polar-only limit, for example in the extensile case, local polarity alignment can stabilize the system only when [2]

$$\kappa k^2 > \begin{cases} \frac{1}{2}(1-\nu) & \text{for } \nu \le 1/2 \text{ and} \\ \frac{1}{(8\nu)} & \text{for } \nu > 1/2. \end{cases} \qquad (47)$$

For the full system with interacting scalar and polar fields, we derive the minimal diffusion term required to stabilize the system in appendix C. For instance, for the gradient-extensile case we find that diffusion can stabilize the system only if

$$\beta + \kappa k^2 > \begin{cases} \frac{1}{2}(1-\nu) & \text{for } \nu \le 1/2 \text{ and} \\ \frac{1}{(8\nu)} & \text{for } \nu > 1/2. \end{cases} \qquad (48)$$

This means that for diffusion to stabilize the system, the gradient alignment $\beta$ needs to be strong enough to dominate the polar field dynamics. If equation (48) is fulfilled, a finite value of the diffusion constant can stabilize the system: $Dk^2 > D_{\text{crit}}k^2 > 0$ (the analytical expression for $D_{\text{crit}}k^2$ is reported in appendix C). Thus, while those systems that fulfill equation (48) can be stabilized by diffusion alone (i.e. with $\kappa = 0$), a system can never be stabilized by polarity alignment alone (i.e. with $D = 0$). This is due to the instability that arises from the coupling of scalar and polar fields (see adiabatic limit discussed in section 5.1.1). We find that this is also true in the gradient-contractile case (appendix C).

In figure 6, we report the magnitude of $D_{\text{crit}}k^2$ for varying $\beta$, $\nu$, and $\kappa k^2$. For instance, for $\nu = 0$ and without polarity alignment, $\kappa = 0$, the critical diffusion $D_{\text{crit}}k^2$ decreases monotonically with the alignment of polarity to the scalar field, $\beta$ (blue solid curve in figure 6(a)), which intuitively makes sense since diffusion and $\beta$ both act together to stabilize the instability of the polar-only system. However, for larger values of the polarity alignment $\kappa k^2$, we find that $D_{\text{crit}}k^2$ can also monotonically increase with $\beta$ (red and green solid curves in figure 6(a)). In this regime, the polar-only limit is not unstable any more, and diffusion is required to compensate for the instability that arises from the coupling of scalar and polar field (section 5.1.1).

Finally, we summarize the stability criteria for the gradient-extensile and gradient-contractile cases in figures 6(b) and (c), respectively. As expected, the critical diffusion term is generally larger in the gradient-contractile case than in the gradient-extensile case.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Stabilization of the system by diffusion and local polarity alignment, for fixed system size with fixed polarity norm.</figcaption></figure>

(a) Minimal diffusion $D_{\text{crit}}k^2$ required to stabilize a gradient-extensile system depending on $\beta$ and polarity alignment $\kappa k^2$. (b),

(c) Minimal diffusion $D_{\text{crit}}k^2$ required to stabilize the system depending on $\beta$ and shear alignment $\nu$ for (b) gradient-extensile and

(c) gradient-contractile systems. In the white regions, the system cannot be stabilized by diffusion. Polarity alignment is set to $\kappa = 0$. The color bar of panel c applies to panel b as well.

## 7. Discussion

The homogeneously deforming state of polar or nematic wet active matter is subject to the well-known Simha-Ramaswamy instability [1, 2]. This raises the question how, despite this instability, active anisotropic tissue deformation can be robust during animal morphogenesis [8]. Animal morphogenesis is known to be organized by large-scale protein concentration patterns (e.g. morphogen gradients) [8]. Under which conditions could such patterns stabilize anisotropic tissue deformation? To address this question, we examined whether a scalar field gradient can stabilize the homogeneously deforming state of an active polar material, even when the scalar field is advected by material flows.

Focusing on the limit of a large system size, we showed that the homogeneously deforming state is always unstable in gradient-contractile systems, i.e. when the active anisotropic stress is, mediated by the polarity, contractile along the direction of the gradient. However, the system can be marginally stable in the gradient-extensile case. This is true both when the active anisotropic stress is controlled directly by the scalar field gradient (section 3, figure 2) and when this control is mediated through the polar field (section 5, figure 4). Intriguingly, in the biology literature we found many instances of animal morphogenetic systems where the effective coupling between controlling morphogen gradient and tissue deformation is gradient-extensile [12, 13, 41–46]. However, so far we could identify almost no instance of a gradient-contractile coupling. This predominance of gradient-extensile systems in developing animals has been remarked once before in the biology literature [47], and our results provide an explanation purely based on active matter physics.

The only example we have found where something akin to a gradient-contractile coupling has been proposed is convergent-extension of the notochord in the ascidian *Ciona intestinalis* [48]. However, evidence is sparse, and the deforming tissue is very small consisting only of 40 cells, which would facilitate other mechanisms of stabilization (see below).

In our study, gradient-extensile systems can only become *marginally* stable. They do not become strictly stable because for perturbations parallel to the gradient direction, the scalar field decouples from flow, and scalar field perturbations along this direction are marginal, i.e. they do neither grow nor shrink. This will be different in more realistic systems, where protein gradients are created, e.g. by secretion, diffusion, and degradation [8, 49, 50].

We further showed that the stability of gradient-extensile systems strongly depends on how the polarity *magnitude* is controlled, where we compare the cases of finite and ZPM (section 5, figure 4, 5). Both kinds have examples in biological tissues. Finite-preferred-magnitude (FPM) polarity resembles Core/Frizzled planar cell polarity (PCP), which is believed to emerge without any cell-external cues [51]. Meanwhile, zero-preferred-magnitude (ZPM) polarity resembles Fat PCP [52] and actin polarity in the *Drosophila* germ band [22, 24], which are believed to show no polarization in the absence of external cues. We demonstrate that for FPM, the response of polarity perturbations to scalar field perturbations is more anisotropic than for ZPM, which affects the generated active stresses and material flows. As a consequence, systems with ZPM are generally more stable than systems with FPM. However, in both cases the coupling of the polar field to the scalar field gradient generally has a stabilizing effect.

We also demonstrated the emergence of parameter regions with transiently growing perturbation modes. These arise because for systems under pure shear deformations the solutions to the linearized dynamics are
co-deforming modes, whose wave vectors are changing during the overall system deformation. As a consequence, the growth rate of the perturbation modes can change over time. We identified parameter regions where some perturbation modes decay, while for other modes the amplitude only transiently grows before decaying forever. To characterize the stability of the system in these regions, we introduce a maximal amplification factor S, which these modes experience while their amplitude grows.

While we discussed here one potential mechanism to stabilize the Simha-Ramaswamy instability during animal morphogenesis, there may be other mechanisms as well. First, in most of our article we focused on the infinite system size limit. For finite polar-only systems, polarity alignment *K*, possibly in combination with the boundary conditions for the polarity, can also stabilize the system. There is a maximal tissue size scale $L_c \sim \sqrt{K/\gamma\tilde{v}}$, where $\tilde{v}$ is the free tissue deformation rate, beyond which polarity alignment becomes insufficient to stabilize the system [2] (section 6). Combining known orders of magnitude for active fly wing deformation $\tilde{v} \sim (10^{-2}\dots10^{-1}) \text{ h}^{-1}$ [17], and a measured PCP alignment rate in the wing, $K/\gamma \sim (1\dots10) \text{ } \mu\text{m}^2 \text{ h}^{-1}$ [21], we obtain $L_c \sim 10 \text{ } \mu\text{m}$, which is on a similar order as values measured in monolayers of different mouse cell lines of $L_c \sim 40 \text{ } \mu\text{m}$ [7]. However, many developing tissues have sizes of $\sim(10\dots100) \text{ } \mu\text{m}$, suggesting that an effective cell polarity alignment is not necessarily sufficient to stabilize the Simha-Ramaswamy instability. When including diffusion and polarity alignment in our calculations for the full system of interacting scalar and polar fields, we find that when scalar and polar fields interact, polarity alignment alone is never sufficient to stabilize the system, at least for fixed polarity norm, while diffusion can be (section 6). This is interesting given that there are biological systems where diffusion is not believed to occur, such as during *Drosophila* germ band extension.

Second, the Simha-Ramaswamy instability occurs only in wet active matter, i.e. when momentum (and angular momentum [53]) is conserved. In other words, active tissue deformation could also be stabilized, e.g. by friction $\zeta$ with a substrate whenever the hydrodynamic length scale $L_h = \sqrt{\eta/\zeta}$ is sufficiently small as compared to $L_c$. This could be relevant for instance for *Drosophila* germ band extension [54]. However, many morphogenetic tissue deformation processes, may be better described as wet active matter. This includes for instance vertebrate limb bud elongation [16] and morphologically similar processes [13], as well as hydra morphogenesis [15].

Third, active oriented materials may be stabilized by lifting the condition of incompressibility [55]. Indeed, developing tissues can show some degree of compressibility. For instance a finite tissue bulk viscosity can arise from cell division and death [56]. However, such a bulk viscosity may become visible only on time scales above the cell division time, which is typically on the order of hours or days. This is also the approximate range of typical anisotropic tissue deformation processes during development. Thus, a bulk viscosity due to cell division would be relevant in particular for slow tissue deformation processes. Moreover, layered 2D tissues called epithelia may additionally exhibit limited 2D compressibility through variation of layer height or cell extrusion. However, a finite compressibility alone is not necessarily sufficient to stabilize deforming systems. For instance, while for a polar-only system with fixed size, a finite bulk viscosity can indeed stabilize the system for sufficiently large value of $\nu$, this stabilizing effect only acts on angles other than multiples of $\pi/2$, i.e. it does not act on bend or splay modes. Consequentially, for deforming systems, a finite bulk viscosity would help stabilize regions with transiently growing modes (e.g. yellow regions in figure 3(j)), but it would not stabilize unstable regions (e.g. red regions in figure 3(j)). For instance, a freely deforming polar-only system is unstable also if the tissue is compressible (green dashed line in figure 3(j)).

Our work prompts for different kinds of experiments on animal morphogenetic systems to test our ideas. First, are really all active anisotropic tissue deformation processes gradient-extensile? So far, there are many systems where the precise role of morphogens for tissue deformation is still unknown [15–18, 57]. In many systems, it is believed that scalar field gradients control the orientation anisotropic deformation, but more experimental evidence is required. Also, is notochord convergent-extension of *Ciona intestinalis* indeed a counter example, i.e. a gradient-contractile system? If so, it is possibly not the only one. How is tissue deformation stabilized in these systems?

Second, relatively little is currently known about what kinds of polarity are used to control tissue deformation, i.e. whether they are ZPM or FPM polarity. In the few systems where more is known, it appears that ZPM polarity controls tissue deformation [12, 41], which in our analysis leads to a more robust behavior. Is this kind of polarity indeed more often used to control tissue deformation in morphogenesis?

Third, many parameter values are still unknown, even in the best-studied biological systems. For instance, while something like a shear alignment effect has been observed in a few systems now [7, 20, 21, 34, 58–61], we are aware of only two systems where the shear alignment parameter $\nu$ has been measured, the *Drosophila* wing [20, 21, 58] and certain cell monolayers *in vitro* [7, 34, 61]. Measuring parameters like this in more developmental systems will allow to quantitatively test our predictions.

While a foundational motivation of active matter physics has always been to better understand collective motion in living systems [62], a lot remains to be learned at the direct interface with biology. Here, we
provide an example for how active matter physics may reveal fundamental principles for animal morphogenesis.

## Data availability statement

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

## Acknowledgments

We thank the Centre Interdisciplinaire de Nanoscience de Marseille (CINaM) for providing office space and S Gsell for a critical reading of the manuscript. The project leading to this publication has received funding from France 2030, the French Government program managed by the French National Research Agency (ANR-16-CONV-0001), and from the Excellence Initiative of Aix-Marseille University—A*MIDEX.

## Appendix A. Co-deforming coordinates

### A.1. Definition

In order to more conveniently solve the linearized dynamics (appendix A.3), we introduce co-deforming coordinates $(\bar{\mathbf{r}}, \bar{\mathbf{t}}) = (\bar{x}, \bar{y}, \bar{t})$. These coordinate map to the lab coordinates $(\mathbf{r}, t) = (x, y, t)$ in the following way:

$$r_i = s_{ij}(\bar{\mathbf{t}})\bar{r}_j \qquad (A1)$$

$$t = \bar{t}, \qquad (A2)$$

where $s(\bar{\mathbf{t}})$ is a time-dependent shear tensor, given by

$$\mathbf{s}(\bar{\mathbf{t}}) = \begin{pmatrix} l_x(\bar{\mathbf{t}}) & 0 \\ 0 & l_x^{-1}(\bar{\mathbf{t}}) \end{pmatrix}. \qquad (A3)$$

Thus, while at some time $t$, lab coordinates range from $0 \le x < L_x(t)$ and $0 \le y < L_y(t)$, co-deforming coordinates map these affinely to the box dimensions at time zero, with $0 \le \bar{x} < L_x(0)$ and $0 \le \bar{y} < L_y(0)$.

As a direct consequence of equations (A1) and (A2), partial derivatives of some quantity $q$ transform as:

$$\bar{\partial}_j q := \frac{\partial q(\bar{\mathbf{r}}, \bar{\mathbf{t}})}{\partial \bar{r}_j} = (\partial_i q) s_{ij} \qquad (A4)$$

$$\bar{\partial}_i q := \frac{\partial q(\bar{\mathbf{r}}, \bar{\mathbf{t}})}{\partial \bar{t}} = \partial_i q + (\partial_i q) \dot{s}_{ij} \bar{r}_j, \qquad (A5)$$

where $\partial_i q := \partial q(\mathbf{r}, t)/\partial r_i$, $\partial_i q := \partial q(\mathbf{r}, t)/\partial t$, and $\dot{s}_{ij} := ds_{ij}/dt = ds_{ij}/d\bar{t}$. Thus, the partial time derivative in co-deforming coordinates, $\partial_i q$, i.e. for fixed $\bar{\mathbf{r}}$, includes a term related to the box shear rate as compared to the partial time derivative with respect to lab coordinates.

Moreover, with the co-deforming Fourier transformation of a quantity $q$ defined as in equation (24), we have the usual derivation rule, where the Fourier transform of $\bar{\partial}_j q(\bar{\mathbf{r}}, \bar{\mathbf{t}})$ is $ik_j q(\bar{k}, \bar{\mathbf{t}})$. From equations (24) and (A1) also follows that a given co-deforming Fourier mode with wave vector $\bar{k}$ corresponds to a lab frame Fourier mode with wave vector $\bar{k}$ with components

$$k_i = \bar{k}_j s_{ji}^{-1}, \qquad (A6)$$

because then we have $\bar{k} \cdot \bar{\mathbf{r}} = \mathbf{k} \cdot \mathbf{r}$.

### A.2. Velocity

To obtain the mapping for the velocity field, we consider a tracer particle that is perfectly advected with the flows. The velocity of that tracer particle corresponds to a total time derivative $v_i = dr_i/dt$, for which we obtain by insertion of equation (A1):

$$v_i = \dot{s}_{ij} \bar{r}_j + s_{ij} \bar{v}_j, \qquad (A7)$$
where $\bar{v}_i = d\bar{r}_i/dt = d\bar{r}_i/d\bar{t}$ is the co-deforming velocity, with $\bar{r}(t)$ being the tracer trajectory in co-deforming coordinates. The first term in equation (A7) corresponds to a motion due to the affine transformation according to box coordinates. Thus, $\bar{v}_i$ can be interpreted as the non-affine component of the flow field.

To obtain a transformation formula for the convective derivative, we consider again our tracer and the presence of some spatio-temporal field $q$. The convective derivative corresponds to the total derivative of the value of $q$ that the tracer locally sees. Thus, we expect analogous expressions for the convective derivative in both lab and co-deforming systems, $\dot{q} := dq/dt = dq/d\bar{t}$. Indeed, using equations (A4)–(A7), we obtain:

$$ \dot{q} = \partial_t q + v_i(\partial_i q) = \bar{\partial}_t q + \bar{v}_i(\bar{\partial}_i q). \qquad (A8) $$

### A.2.1. Dynamical equations

The dynamical equations for scalar and polar fields, equations (12) and (13), in co-deforming coordinates are:

$$ \frac{dc}{d\bar{t}} = 0 \qquad (A9) $$

$$ \frac{dp_i}{d\bar{t}} = -\frac{g(p)}{\tau} p_i - \nu \tilde{v}_{ij}^0 p_j + \beta s_{ij}^{-1} \bar{\partial}_j c - \frac{1}{2} \left[ (\nu + 1) s_{li}^{-1} s_{jk} + (\nu - 1) s_{lj}^{-1} s_{ik} \right] \bar{\partial}_l \bar{v}_k p_j. \qquad (A10) $$

Here, we left out diffusion and polarity alignment, and we use the box shear rate $\tilde{v}_{ij}^0 = (\dot{s}_{jk}s_{ki}^{-1} + \dot{s}_{ik}s_{kj}^{-1})/2$, which is for the box shear tensor defined in equation (A3):

$$ \tilde{v}^0 = \begin{pmatrix} \dot{l}_x/l_x & 0 \\ 0 & -\dot{l}_x/l_x \end{pmatrix}. \qquad (A11) $$

We do not rewrite in co-deforming coordinates the incompressible Stokes' equations, equations (14) and (15), because this will not be required in what follows.

### A.3. Linearized dynamics around the homogeneously deforming state

We linearize the dynamics around the homogeneously deforming state, given by

$$ c_0 = \bar{x} \qquad p_0 = p_0 \hat{\mathbf{x}} \qquad \bar{v}_0 = \mathbf{0}. \qquad (A12) $$

Using co-deforming coordinates will facilitate dealing with the advective terms when solving the linearized dynamics.

To fix a value for $p_0$ at some time point $t$, we use equations (A10) and (A12) with the stationarity condition $dp_x/dt = 0$:

$$ g(p_0) = \frac{\beta \tau}{p_0 l_x(t)} - \nu \tilde{v}_{xx}^0 \tau. \qquad (A13) $$

For constant $\beta$ and a deforming box, the state $p_0$ given by this equation is only transiently stationary, due to the time-dependent $l_x$. However, in the adiabatic limit where polarity relaxation is much faster than box deformation, $\tau |\tilde{v}_{xx}^0| \ll 1$, the homogeneous dynamics will generally be close to the state $p_0$ given by equation (A13). In the main text, we circumvent these issues by setting

$$ \beta(t) = \beta_0 l_x(t). \qquad (A14) $$

with constant $\beta_0$. In this case equation (A13) always defines a stationary solution of the dynamics.

We consider the following perturbation to linear order:

$$ c = \bar{x} + \delta c \qquad (A15) $$

$$ \mathbf{p} = p_0 \hat{\mathbf{x}} + \delta \mathbf{p} \qquad (A16) $$

$$ \bar{v} = \delta \bar{v}. \qquad (A17) $$

Insertion into equations (A9) and (A10) yields to first order in co-deforming Fourier space:

$$ \bar{\partial}_t \delta c(\bar{k}, \bar{t}) = -\delta \bar{v}_x \qquad (A18) $$
$$ \bar{\partial}_t \delta p_x(\bar{k}, \bar{t}) = \left( -\frac{1}{\tau} \left[ g(p_0) + g'(p_0)p_0 \right] - \nu \tilde{v}_{xx}^0 \right) \delta p_x + i\beta l_x^{-1} \bar{k}_x \delta c - i\nu p_0 \bar{k}_x \delta \bar{v}_x \quad (A19) $$

$$ \bar{\partial}_t \delta p_y(\bar{k}, \bar{t}) = \left( -\frac{1}{\tau} g(p_0) + \nu \tilde{v}_{xx}^0 \right) \delta p_y + i \beta l_x \bar{k}_y \delta c - \frac{i p_0}{2} \left[ (\nu + 1) l_x^2 \bar{k}_y \delta \bar{v}_x + (\nu - 1) l_x^{-2} \bar{k}_x \delta \bar{v}_y \right]. \quad (A20) $$

Using the co-deforming coordinates allowed us to simplify the advective terms, which otherwise lead to spatially dependent coefficients. Here, for simplicity, we left out the diffusion and local polarity alignment terms.

To close the system of equations, equations (A18)–(A20), we need to obtain an expression for $\delta \bar{v}$ by solving the incompressible Stokes’ equations, equations (14) and (15). We insert the expression for the lab-frame velocity equation (A7) into these equations, and obtain after lab-frame Fourier transformation of the linearized dynamics:

$$ \delta \bar{v}_m = \frac{is_{mi}^{-1} k_l}{\eta k^2} \left( \delta_{ij} - \frac{k_i k_j}{k^2} \right) \delta \tilde{\sigma}_{lj}^a. \quad (A21) $$

Here, $k$ is the magnitude of the lab-frame wave vector $k$. Using that to linear order the active stress $\delta \tilde{\sigma}_{lj}^a$ is given by equation (30), we obtain:

$$ \delta \bar{v}_x = \frac{il_x^{-1}}{k} \operatorname{sgn}(\alpha) p_0 \sin(\phi) \delta \Phi \quad (A22) $$

$$ \delta \bar{v}_y = -\frac{il_x}{k} \operatorname{sgn}(\alpha) p_0 \cos(\phi) \delta \Phi \quad (A23) $$

with $\delta \Phi = \sin(2\phi)\delta p_x - \cos(2\phi)\delta p_y$, and $\phi$ being the angle of the lab-frame wave vector $k$. Note that a transformation of equations (14) and (15) into co-deforming coordinates, followed by a linearization and a co-deforming Fourier transformation leads to the same result.

Inserting the velocity perturbation, equations (A22) and (A23), into the linearized dynamics, equations (A18)–(A20), we obtain:

$$ \bar{\partial}_t \delta c = -Dk^2 \delta c - \frac{il_x^{-1}}{k} \operatorname{sgn}(\alpha) p_0 \sin \phi \sin 2\phi \delta p_x + \frac{il_x^{-1}}{k} \operatorname{sgn}(\alpha) p_0 \sin \phi \cos 2\phi \delta p_y \quad (A24) $$

$$ \bar{\partial}_t \delta p_x = i\beta k \cos \phi \delta c + \left[ -\frac{1}{\tau} \left( g(p_0) + g'(p_0)p_0 \right) - \nu \tilde{v}_{xx}^0 + \frac{\operatorname{sgn}(\alpha) p_0^2}{2} \nu \sin^2 2\phi - \kappa k^2 \right] \delta p_x \\ - \frac{\operatorname{sgn}(\alpha) p_0^2}{2} \nu \sin 2\phi \cos 2\phi \delta p_y \quad (A25) $$

$$ \bar{\partial}_t \delta p_y = i\beta k \sin \phi \delta c - \frac{\operatorname{sgn}(\alpha) p_0^2}{2} \sin 2\phi (\nu \cos 2\phi - 1) \delta p_x \\ + \left[ -\frac{1}{\tau} g(p_0) + \nu \tilde{v}_{xx}^0 + \frac{\operatorname{sgn}(\alpha) p_0^2}{2} \cos 2\phi (\nu \cos 2\phi - 1) - \kappa k^2 \right] \delta p_y. \quad (A26) $$

The derivatives on the left-hand sides are partial derivatives for constant co-deforming wave vectors $\bar{k}$. Thus, solutions to the linearized dynamics are co-deforming Fourier modes with time-dependent amplitude. The right-hand side is written in terms of angle $\phi$ and magnitude $k$ of the lab-frame wave vector out of convenience only.

## Appendix B. Linear stability of the homogeneously deforming state

Here, we discuss the linear stability of the homogeneously deforming state based on the linearized dynamics, equations (A24)–(A26).
### B.1. Gradient-contractile systems ($\operatorname{sgn}(\alpha) = 1$) are always unstable

Here we show that the system is always unstable for $\operatorname{sgn}(\alpha) = 1$ and $\beta > 0$, and for vanishing diffusion and polarity alignment. This is because in this case there is a co-deforming mode with wave vector direction $\tilde{\phi} = \pi/2$ that is permanently growing. For this specific angle, this corresponds to a lab-frame wave vector with the same angle, $\phi = \pi/2$.

To prove that there is a mode with $\phi = \pi/2$ that is permanently growing, we show that for this angle the characteristic polynomial $P(\omega)$ of the matrix describing the linearized dynamics, equations (A24)–(A26), has at least one positive zero. For $\phi = \pi/2$ this polynomial is:

$$ P(\omega) = \operatorname{sgn}(\alpha)\beta l_x^{-1} \left[ \beta l_x^{-1} + \frac{1}{\tau}g'(p_0)p_0^2 \right] + \operatorname{sgn}(\alpha)\beta l_x^{-1}p_0\omega - \omega \left[ -\frac{1}{\tau}\left(g(p_0) + g'(p_0)p_0\right) - \nu\tilde{v}_{xx}^0 - \omega \right] \times \left[ -\frac{1}{\tau}g(p_0) + \nu\tilde{v}_{xx}^0 + \frac{\operatorname{sgn}(\alpha)p_0^2}{2}(\nu+1) - \omega \right], \quad (B1) $$

where we also used equation (A13) to simplify the absolute-order term in $\omega$ (first term on the right-hand side).

The polynomial $P(\omega)$ has at least one positive root, because first, $P(0) > 0$, since the absolute term is positive for $\operatorname{sgn}(\alpha) = 1$. Second, the coefficient in front of the cubic term in $\omega$ is negative, so that $P(\omega) \to -\infty$ for $\omega \to \infty$. Thus, using the intermediate value theorem, it follows that $P(\omega)$ has at least one positive zero. As a consequence, any gradient-contractile system is always unstable.

Of course, this does not preclude that it could in principle be possible to stabilize gradient-contractile systems when including diffusion and/or polarity alignment in a system with a finite size.

### B.2. Scalar field only

Here we briefly discuss the limit where $\operatorname{sgn}(A) = 1$, $B=0$, i.e. $g(p) = 1$, and $\tau \to 0$ while $\beta\tau = 1$. In this limit, polarity relaxation is much faster than box deformation, and equation (A13) is solved by

$$ p_0 = l_x^{-1}. \quad (B2) $$

Moreover, polarity perturbation away from this state relaxes adiabatically fast towards the scalar field perturbation. Using the linearized dynamics for the polarity, equations (A25) and (A26):

$$ \delta p_x = ik \cos \phi \delta c \quad (B3) $$

$$ \delta p_y = ik \sin \phi \delta c. \quad (B4) $$

Insertion in the linearized dynamics of the scalar field, equations (A24), yields:

$$ \bar{\partial}_t \delta c = \operatorname{sgn}(\alpha) l_x^{-2} \sin^2 \phi \delta c. \quad (B5) $$

### B.3. Polar field only

Here, we discuss the case where $\beta = 0$, i.e. where the scalar field is irrelevant. We focus on the extensile case, $\operatorname{sgn}(\alpha) = -1$.

#### B.3.1. Fixed magnitude, $\operatorname{sgn}(A) = -1, \tau \to 0$

In the case of fixed polarity magnitude, $\operatorname{sgn}(A) = -1$ and $\tau \to 0$, we have $p_0 = 1$ and $\delta p_x = 0$. In this limit, equations (A25) and (A26) yield:

$$ \bar{\partial}_t \delta p_y = \omega(\phi) \delta p_y \quad \text{with} \quad (B6) $$

$$ \omega(\phi) = 2\nu\tilde{v}_{xx}^0 - \frac{1}{2} \left[ \nu \cos^2(2\phi) - \cos 2\phi \right]. \quad (B7) $$

In order to plot the analytical boundaries of the phase diagram in figure 3(j), we study the growth rates of (a) the bend mode, $\omega(\phi = 0)$, (b) the splay mode $\omega(\phi = \pi/2)$, and (c) the maximum growth rate $\omega$ maximized over all angles, $\omega_{max}$. When at least one of the bend or the splay mode has a positive growth rate, the system is unstable. Otherwise, when both modes have negative growth rates, the system is either stable when $\omega_{max}$ is negative, or stable with transiently growing modes whenever $\omega_{max}$ is positive.
<figure><img src="image_10.png" /><figcaption>Figure B7. Maximal perturbation growth rate, $\omega_{max}$, for a polar-only system with fixed system size and finite preferred polarity magnitude. We plot three different values of polarity relaxation timescale, $\tau$. As $\tau \to +\infty$, the growth rates in the $\nu > 1$ regime converge to zero.</figcaption></figure>

The growth rates of bend and splay modes are:

$$ \omega(\phi = 0) = 2\nu\tilde{v}_{xx}^0 - \frac{1}{2}(\nu - 1) \qquad (B8) $$

$$ \omega(\phi = \pi/2) = 2\nu\tilde{v}_{xx}^0 - \frac{1}{2}(\nu + 1). \qquad (B9) $$

The maximum of $\omega(\phi)$ over all angles $\phi$ is:

$$ \omega_{\max} = \begin{cases} \omega(\phi = 0) & \text{for } \nu \le 1/2 \\ 2\nu\tilde{v}_{xx}^0 + 1/(8\nu) & \text{for } \nu > 1/2. \end{cases} \qquad (B10) $$

In the phase diagram in figure 3(j), there are two kinds of solid curves which define the boundaries of regions of different behavior. The first kind of curve, defined over the whole $\nu$ range, satisfies the $\omega(\phi = 0) = 0$ equation. This curve defines the boundaries of unstable regions as $\omega(\pi/2) < \omega(0)$ for every $(\nu, \tilde{v}_{xx}^0)$. The second kind of curve, only present for $\nu > 0.5$, is the $\omega_{max} = 0$ curve, which defines the region of transiently growing modes. Above this curve, both bend and splay modes have negative growth rates, but the maximal growth rate is positive.

### B.3.2. Zero preferred magnitude, sgn(A) = 1

Here we derive analytical expressions for the phase boundaries in figure 3(m). For zero preferred polarity magnitude, a finite polarity magnitude

$$ p_0 = \sqrt{-1 - \nu\tilde{v}_{xx}^0\tau} \qquad (B11) $$

exists only for large enough box shear rate $|\tilde{v}_{xx}^0| > 0$, in particular $\nu\tilde{v}_{xx}^0\tau < 1$. When this condition is not met, the polarity in the stationary state has zero magnitude.

Diagonalizing the $2 \times 2$ matrix that corresponds to the dynamics of $\delta p_x$ and $\delta p_y$, equations (A25) and (A26), we find the $\phi$-dependent amplitude growth rates $\omega$:

$$ \omega_{1,2} = \frac{1}{2} \left[ \mathcal{T} \pm \sqrt{\mathcal{T}^2 - 4\mathcal{D}} \right], \quad \text{where} \qquad (B12) $$

$$ \mathcal{T}(\phi) = \frac{p_0^2}{2}(\cos 2\phi - \nu) - \frac{2}{\tau}(1 + 2p_0^2) \qquad (B13) $$

$$ \mathcal{D}(\phi) = \frac{p_0^2}{\tau} \left[ \left( \nu + \frac{4}{\tau} \right) (1 + p_0^2) - \cos 2\phi \left( \nu \cos 2\phi + p_0^2 \right) \right]. \qquad (B14) $$

Based on equation (B12), a mode with angle $\phi$ grows in amplitude only if $\mathcal{T}(\phi) > 0$ or $\mathcal{D}(\phi) < 0$.
As we argued in the section before, a system is in the unstable regime whenever the bend ($\phi = 0$) or the splay ($\phi = \pi/2$) mode grows. To see when this is the case, we first remark that $T(0) > T(\pi/2)$ while $D(0) < D(\pi/2)$. Thus, the bend mode will always be growing when the splay mode is, and so the bend mode is sufficient to decide whether the system is in the unstable regime. Moreover, we find that $D(0) = p_0^2[-2T(0) - 4p_0^2/\tau]/\tau$, from which follows that whenever $D(0) \ge 0$ then $T(0) \le 0$. Taken together, this means that the system is in the unstable regime if and only if $D(0) < 0$. Indeed, we find that the criterion $D(0) = 0$ defines the boundaries of all unstable regimes in figure 3(m), where for a given parameter point ($\nu, \tilde{v}_{xx}^0$), the polarity magnitude $p_0$ needs to be inserted from equation (B11).

Further, if the system is not unstable, for $D(0) \ge 0$, the regime could be either stable or stable with transiently growing modes. To know if there are transiently growing modes, we need to know if there is any $\phi$ with a positive $\omega$, i.e. with $T(\phi) > 0$ or $D(\phi) < 0$. However, we have $D(0) \ge 0$ and thus $T(0) \le 0$. Moreover, the maximum of $T(\phi)$ is at $\phi = 0$. Thus, we have $T(\phi) \le 0$ for any $\phi$ outside the unstable regime. As a consequence, any regime with transiently growing modes needs to have an angle $\phi$ for which $D(\phi) < 0$. Whether this is the case depends on the sign of $\nu$. First, for $\nu \ge 0$, the minimum of $D(\phi)$ is at $\phi = 0$. Thus, outside of the unstable regime, where $T(0) < 0$ and thus $D(0) > 0$, the value of $D(\phi)$ can never be negative for any $\phi$. Thus, for $\nu \ge 0$, there are no stable regimes with transiently growing modes (see figure 3(m)). Second, for $\nu < 0$, it can be shown that the minimum of $D(\phi)$ is not at $\phi = 0$ only if $p_0^2 < -2\nu$. In these cases, the minimal value of $D(\phi)$ is

$$D_{\min} = \frac{p_0^2}{\tau} \left[ \left( \nu + \frac{4}{\tau} \right) (1 + p_0^2) + \frac{p_0^4}{2\nu} \right]. \qquad (B15)$$

Taken together, for $T(0) < 0$ and $\nu < 0$ the system has transiently growing modes only if $p_0^2 < -2\nu$ and $D_{\min} < 0$.

## B.4. Scalar and polar field

In this section we focus exclusively on the extensile case, $\operatorname{sgn}(\alpha) = -1$.

### B.4.1. Fixed system size, strong coupling to the scalar field

Here we discuss the limit of strong coupling to the scalar field. To more clearly understand what are the correct parameter regimes for strong and weak coupling to the scalar field gradient, we re-express the linearized dynamics using the following definitions:

$$\delta \tilde{c} := ikp_0 \delta c \qquad (B16)$$

$$\tilde{\beta} := \frac{\beta_0}{p_0^3} = \frac{g(p_0)}{\tau p_0^2}. \qquad (B17)$$

The last equation on the second line follows directly from the stationary-state requirement for $p_0$, equation (A13), for fixed system size. We now use these two definitions to simplify the linearized dynamics, equations (A24)–(A26):

$$\frac{1}{p_0^2} \bar{\partial}_t \delta \tilde{c} = - \sin \phi \sin 2\phi \delta p_x + \sin \phi \cos 2\phi \delta p_y \qquad (B18)$$

$$\frac{1}{p_0^2} \bar{\partial}_t \delta p_x = \tilde{\beta} \cos \phi \delta \tilde{c} - \left[ \tilde{\beta} G^{-1} + \frac{\nu}{2} \sin^2 2\phi \right] \delta p_x + \frac{\nu}{2} \sin 2\phi \cos 2\phi \delta p_y \qquad (B19)$$

$$\frac{1}{p_0^2} \bar{\partial}_t \delta p_y = \tilde{\beta} \sin \phi \delta \tilde{c} + \frac{1}{2} \sin 2\phi (\nu \cos 2\phi - 1) \delta p_x - \left[ \tilde{\beta} + \frac{1}{2} \cos 2\phi (\nu \cos 2\phi - 1) \right] \delta p_y. \qquad (B20)$$

In equation (B19), we have also used the definition of $G(p_0)$, equation (43). From equations (B19) and (B20), we see directly that $\tilde{\beta}$ compares the scalar field coupling strength (and the polarity magnitude control) with the flow-induced feedback. Thus, the correct limit for a strong scalar field coupling is $\tilde{\beta} \gg 1$, i.e. $\beta \gg p_0^3$.

For strong scalar field coupling, $\tilde{\beta} \gg 1$, equations (B19) and (B20) imply that polarity relaxes adiabatically fast towards:

$$\delta p_x = iGp_0 k_x \delta c \qquad (B21)$$

$$\delta p_y = ip_0 k_y \delta c. \qquad (B22)$$
These equations can be used to obtain a criterion for the stability of the system. Indeed, combining equations (A18), (A21), (30), (B21) and (B22), we obtain:

$$ \partial_t \delta c = p_0^2 \sin^2 \phi [2(1-G)\cos^2 \phi - 1] \delta c. \qquad (B23) $$

From this equation directly follows that the system is marginally stable whenever $G \ge 1/2$.

Alternatively, to allow for the more intuitive explanation in the main text, we introduce an angle $\theta_p$, which globally characterizes the orientation of $\delta p$. However, for a given perturbation mode with wave vector $k$, the direction of $\delta p$ will spatially depend on the phase of the Fourier mode. We remove this ambiguity by dividing $\delta p$ by $ik\delta c$, defining the angle $\theta_p$ by

$$ \frac{\delta \mathbf{p}}{ik\delta c} = \hat{p} \left( \frac{\cos \theta_p}{\sin \theta_p} \right), \qquad (B24) $$

where $\hat{p} > 0$. From equations (B21), (B22) and (B24) then follows for the norm $\hat{p} = p_0 \hat{G}$, where $\hat{G} = [G^2 \cos^2 \phi + \sin^2 \phi]^{1/2}$, and for the angle:

$$ \tan \theta_p = \frac{1}{G} \tan \phi. \qquad (B25) $$

This is equation (44) in the main text. We proceed similarly to define the angle of the active stress perturbation nematic, $\theta_\sigma$. More precisely, we define $\theta_\sigma$ as the angle of the nematic $-\delta\tilde{\sigma}_{ij}/(ik\delta c)$:

$$ -\frac{\delta \tilde{\sigma}^a}{ik\delta c} = \hat{\sigma} \begin{pmatrix} \cos 2\theta_\sigma & \sin 2\theta_\sigma \\ \sin 2\theta_\sigma & -\cos 2\theta_\sigma \end{pmatrix}, \qquad (B26) $$

where $\hat{\sigma} > 0$. Together with equations (30) and (B24), we have indeed

$$ \theta_\sigma = \frac{\theta_p}{2}. \qquad (B27) $$

and $\hat{\sigma} = -\alpha p_0 \hat{p} = -\alpha p_0^2 \hat{G}$. Insertion of equation (B26) into the equation for the velocity perturbation $\delta \bar{v}_x$, equation (A21), and noting that for fixed system size $\delta \bar{v}_x = \delta v_x$, yields:

$$ \delta v_x = -p_0^2 \hat{G} \sin \phi \sin (2[\theta_\sigma - \phi]) \delta c. \qquad (B28) $$

Together with equation (A18), this results in

$$ \partial_t \delta c = p_0^2 \hat{G} \sin \phi \sin (2[\theta_\sigma - \phi]) \delta c. \qquad (B29) $$

Hence, the system is unstable whenever there is a $\phi$ with positive $\sin \phi \sin (2[\theta_\sigma - \phi])$.

### B.4.2. Deforming system

In this part we derive the analytical curves that define the unstable regions in figure 5. As discussed in appendix section B.3, the system is in the unstable regime, whenever a mode grow with $\phi = 0$ or with $\phi = \pi/2$ grows.

For $\phi = 0$, the linearized dynamics of $\delta p_y$, (A26), decouple from those of $\delta c$ and $\delta p_x$, equations (A24) and (A25). We find for the growth rate of the orientational perturbations, $\delta p_y$:

$$ \omega(\phi = 0) = -\frac{\beta_0}{p_0^2} + 2\nu \tilde{v}_{xx}^0 - \frac{p_0^2}{2}(\nu - 1). \qquad (B30) $$

Moreover, for $\phi = \pi/2$, the linearized dynamics of $\delta c$ and $\delta p_y$ decouples from that of $\delta p_x$, where it can be directly shown that the maximum growth rate of the $\phi = \pi/2$ modes is always smaller than $\omega(\phi = 0)$. Thus, the $\phi = 0$ growth rate decides whether the system is in the unstable regime or not. Hence, all boundaries to the unstable regime are given by $\omega(\phi = 0) = 0$ (all black solid curves in figure 5). To obtain a relation between $\nu$ and $\tilde{v}_{xx}^0$, for the fixed norm case, $p_0$ was set to one (figures 5(a)–(c)), while for ZPM polarity, the polarity magnitude $p_0$ was eliminated using equation (40) (figures 5(d)–(f)).
## Appendix C. System size effects

Here we report analytical results for the minimal diffusion term, $D_{crit}k^2$ required to stabilize the system, depending on the parameters $\beta$ and $\nu$, the magnitude of the polarity alignment term $\kappa k^2$, and the sign of the activity sgn($\alpha$). We focus on a non-deforming system with fixed polarity magnitude. As discussed in the main text, without diffusion or polarity alignment the system is always unstable in this case.

The linearized perturbation dynamics read

$$ \bar{\partial}_t \delta c = -Dk^2 \delta c + \frac{\operatorname{isgn}(\alpha)}{k} \sin\phi \cos 2\phi \delta p_y \quad (C1) $$

$$ \bar{\partial}_t \delta p_y = i\beta k \sin\phi \delta c + \left[ \frac{\operatorname{sgn}(\alpha)}{2} \cos 2\phi (\nu \cos 2\phi - 1) - \beta - \kappa k^2 \right] \delta p_y. \quad (C2) $$

The perturbation growth rates $\omega_{1,2}(\phi)$ are the eigen values of the coefficient matrix. Denoting trace and determinant of this matrix by $\mathcal{T}(\phi)$ and $\mathcal{U}(\phi)$, respectively, we have:

$$ \omega_{1,2}(\phi) = \frac{1}{2} \left[ \mathcal{T}(\phi) \pm \sqrt{\mathcal{T}^2(\phi) - 4\mathcal{U}(\phi)} \right]. \quad (C3) $$

The system is stable only iff both $\omega_1$ and $\omega_2$ are negative, i.e. iff $\mathcal{T} < 0$ and $\mathcal{U} > 0$, for all angles $\phi$.

As a consequence, we find after some calculation that for gradient-extensile activity, sgn($\alpha$) = -1, the system can be stabilized by diffusion iff

$$ \beta + \kappa k^2 > \begin{cases} (1 - \nu)/2 & \text{for } \nu \le 1/2 \text{ and} \\ 1/8\nu & \text{for } \nu > 1/2. \end{cases} \quad (C4) $$

In both cases, the critical value of the diffusion term is given by

$$ D_{\text{crit}} k^2 = \beta \left[ -1 + 4(\beta + \kappa k^2) + 4 \sqrt{(\beta + \kappa k^2)(\beta + \kappa k^2 - [1 - \nu]/2)} \right]^{-1}. \quad (C5) $$

For gradient-contractile activity, sgn($\alpha$) = 1, the system can be stabilized by diffusion iff

$$ \beta + \kappa k^2 > \begin{cases} (1 + \nu)/2 & \text{for } \nu \ge -1/2 \text{ and} \\ -1/8\nu & \text{for } \nu < -1/2. \end{cases} \quad (C6) $$

For $\nu \ge -1/2$ or whenever $\nu < -1/2$ and $\beta + \kappa k^2 > (1 - \nu)/6$ the critical value of the diffusion term is

$$ D_{\text{crit}} k^2 = \frac{\beta}{\beta + \kappa k^2 - [1 + \nu]/2}. \quad (C7) $$

Meanwhile, whenever $\nu < -1/2$ and $(1 - \nu)/6 > \beta + \kappa k^2 > -1/8\nu$ the critical value of the diffusion term is

$$ D_{\text{crit}} k^2 = \beta \left[ -1 - 4(\beta + \kappa k^2) + 4 \sqrt{(\beta + \kappa k^2)(\beta + \kappa k^2 + [1 - \nu]/2)} \right]^{-1}. \quad (C8) $$

## ORCID iDs

Muhamet Ibrahimi https://orcid.org/0000-0002-8117-8614

Matthias Merkel https://orcid.org/0000-0001-9118-1270

## References

[1] Simha R A and Ramaswamy S 2002 Hydrodynamic fluctuations and instabilities in ordered suspensions of selfpropelled particles Phys. Rev. Lett. 89 058101

[2] Voituriez R, Joanny J F and Prost J 2005 Spontaneous flow transition in active polar gels Europhys. Lett. 70 404

[3] Voituriez R, Joanny J F and Prost J 2006 Generic phase diagram of active polar films Phys. Rev. Lett. 96 1

[4] Sanchez T, Chen D T, Decamp S J, Heymann M and Dogic Z 2012 Spontaneous motion in hierarchically assembled active matter Nature 491 431

[5] Be'er A and Ariel G 2019 A statistical physics view of swarming bacteria Mov. Ecol. 7 1

[6] Liu S, Shankar S, Marchetti M C and Wu Y 2021 Viscoelastic control of spatiotemporal order in bacterial active matter Nature 590 80
[7] Duclos G, Blanch-Mercader C, Yashunsky V, Salbreux G, Joanny J-F, Prost J and Silberzan P 2018 Spontaneous shear flow in confined cellular nematics *Nat. Phys.* **14** 728

[8] Wolpert L, Tickle C and Arias A M 2015 *Principles of Development* (Oxford: Oxford University Press)

[9] Tada M and Heisenberg C-P 2012 Convergent extension: using collective cell migration and cell intercalation to shape embryos *Development* **139** 3897

[10] Sutherland A, Keller R and Lesko A 2019 Convergent extension in mammalian morphogenesis *Semin. Cell Dev. Biol.* **100** 199–211

[11] Bertet C, Sulak L and Lecuit T 2004 Myosin-dependent junction remodelling controls planar cell intercalation and axis elongation *Nature* **429** 667

[12] Bosveld F *et al* 2012 Mechanical control of morphogenesis by Fat/Dachsous/Four-jointed planar cell polarity pathway *Science* **336** 724

[13] Tao H *et al* 2019 Oscillatory cortical forces promote three dimensional cell intercalations that shape the murine mandibular arch *Nat. Commun.* **10** 1703

[14] Keller R and Sutherland A 2020 *Current Topics in Developmental Biology* vol 136 1st edn (Amsterdam: Elsevier Inc.) pp 271–317

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

[16] Sermeus Y, Vangheel J, Geris L, Smees B and Tylzanowski P 2022 Mechanical regulation of limb bud formation *Cells* **11** 1

[17] Etournay R *et al* 2015 Interplay of cell dynamics and epithelial tension during morphogenesis of the Drosophila pupal wing *eLife* **4** e07090

[18] Dye N A, Popovic M, Venkatesan Iyer K, Fuhrmann J F, Piscitello-Gómez R, Eaton S and Julicher F 2021 Self-organized patterning of cell morphology via mechanosensitive feedback *eLife* **10** 1

[19] Zallen J A 2007 Planar polarity and tissue morphogenesis *Cell* **129** 1051

[20] Aigouy B, Farhadifar R, Staple D B, Sagner A, Röper J-C, Jülicher F and Eaton S 2010 Cell flow reorients the axis of planar polarity in the wing epithelium of Drosophil *Cell* **142** 773

[21] Merkel M, Sagner A, Gruber F S, Etournay R, Blasse C, Myers E, Eaton S and Jülicher F 2014 The balance of prickle/spiny-legs isoforms controls the amount of coupling between core and fat PCP systems *Curr. Biol.* **24** 2111

[22] Zallen J A and Wieschaus E 2004 Patterned gene expression directs bipolar planar polarity in Drosophila *Dev. Cell* **6** 343

[23] Butler M T and Wallingford J B 2017 Planar cell polarity in development and disease *Nat. Rev. Mol. Cell Biol.* **18** 375

[24] Lavalou J, Mao Q, Harmansa S, Kerridge S, Lellouch A C, Philippe J M, Audebert S, Camoin L and Lecuit T 2021 Formation of polarized contractile interfaces by self-organized Toll-8/Cirl GPCR asymmetry *Dev. Cell* **56** 1574

[25] Marchetti M C, Joanny J F, Ramaswamy S, Liverpool T B, Prost J, Rao M and Simha R A 2013 Hydrodynamics of soft active matter *Rev. Mod. Phys.* **85** 1143

[26] Jülicher F, Grill S W and Salbreux G 2018 Hydrodynamic theory of active matter *Rep. Prog. Phys.* **81** 076601

[27] Wartlick O, Mumcu P, Kicheva A, Bittig T, Seum C, Jülicher F and González-Gaitán M 2011 Dynamics of Dpp signaling and proliferation control *Science* **331** 1154

[28] Banerjee S, Utuje K J and Marchetti M C 2015 Propagating stress waves during Epithelial expansion *Phys. Rev. Lett.* **114** 1

[29] Blanch-Mercader C, Vincent R, Bazellières E, Serra-Picamal X, Trepat X and Casademunt J 2017 Effective viscosity and dynamics of spreading epithelia: a solvable model *Soft Matter* **13** 1235

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

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

[32] Streichan S J, Lefebvre M, Noll N, Wieschaus E F and Shraiman B I 2018 Global morphogenetic flow is accurately predicted by the spatial distribution of myosin motors *eLife* **7** e27454

[33] Münster S, Jain A, Mietke A, Pavlopoulos A, Grill S W and Tomancak P 2019 Attachment of the blastoderm to the vitelline envelope affects gastrulation of insects *Nature* **568** 395

[34] Blanch-Mercader C, Guillamat P, Roux A and Kruse K 2021 Integer topological defects of cell monolayers: mechanics and flows *Phys. Rev. E* **103** 012405

[35] de Gennes P G and Prost J 1995 *The Physics of Liquid Crystals* (Oxford: Clarendon)

[36] Kruse K, Joanny J F, Jülicher F, Prost J and Sekimoto K 2005 Generic theory of active polar gels: a paradigm for cytoskeletal dynamics *Eur. Phys. J. E* **16** 5

[37] Giomi L, Bowick M J, Mishra P, Sknepnek R and Cristina Marchetti M 2014 Defect dynamics in active nematics *Phil. Trans. R. Soc. A Math. Phys. Eng. Sci.* **372** 20130365

[38] Tiribocchi A, Wittkowski R, Marenduzzo D and Cates M E 2015 Scalar active matter in a momentum-conserving fluid *Phys. Rev. Lett.* **115** 188302

[39] Kirkpatrick T R and Bhattacharyjee J K 2019 Driven active matter: Fluctuations and a hydrodynamic instability *Phys. Rev. Fluids* **4** 1

[40] Giordano M G, Bonelli F, Carenza L N, Gonnella G and Negro G 2021 Activity induced isotropic-polar transition in active liquid crystals *Europhys. Lett.* **133** 58004

[41] Kong D, Wolf F and Großhans J 2017 Forces directing germ-band extension in Drosophila embryos *Mech. Dev.* **144** 11

[42] Johansen K A, Iwaki D D and Lengyel J A 2003 Localized JAK/STAT signaling is required for oriented cell rearrangement in a tubular epithelium *Development* **130** 135

[43] Ninomiya H, Elinson R P and Winklbauer R 2004 Anteroposterior tissue polarity links mesoderm convergent ex tension to axial patterning *Nature* **430** 364

[44] Saxena A, Denholm B, Bunt S, Bischoff M, VijayRaghavan K and Skaer H 2014 Epidermal growth factor signalling controls Myosin II planar polarity to orchestrate convergent extension movements during Drosophila tubulogenesis *PLoS Biol.* **12** 1002013

[45] Alégot H, Pouchin P, Bardot O and Mirouse V 2018 Jak-stat pathway induces Drosophila follicle elongation by a gradient of apical contractility *eLife* **7** e32943

[46] Bénazéraf B, Francois P, Baker R E, Denans N, Little C D and Pourquié O 2010 A random cell motility gradient downstream of FGF controls elongation of an amniote embryo *Nature* **466** 248

[47] Keller R 2006 Mechanisms of Elongation in Embryogenesis *J. Dev.* **133** 2291–302

[48] Shi W, Peyrot S M, Munro E and Levine M 2009 FGF3 in the floor plate directs notochord convergent extension in the Ciona tadpole *Development* **136** 23

[49] Turing A M 1952 The chemical basis of morphogenesis *Phil. Trans. R. Soc. B* **237** 37
[50] Romanova-Michaelides M, Hadjivasiliou Z, Aguilar-Hidalgo D, Basagiannis D, Seum C, Dubois M, Jülicher F and Gonzalez-Gaitan M 2022 *Nature* 602 287–93

[51] Devenport D 2014 Cell biology in development: The cell biology of planar cell polarity *J. Cell Biol.* **207** 171

[52] Brittle A, Thomas C and Strutt D 2012 Planar polarity specification through asymmetric subcellular localization of fat and dachsous *Curr. Biol.* **22** 907

[53] Maitra A, Srivastava P, Marchetti M C, Lintuvuori J S, Ramaswamy S and Lenz M 2018 A nonequilibrium force can stabilize 2D active nematics *Proc. Natl Acad. Sci.* **115** 6934

[54] D'Angelo A, Dierkes K, Carolis C, Salbreux G and Solon J 2019 *In vivo* force application reveals a fast tissue softening and external friction increase during early embryogenesis *Curr. Biol.* **29** 1564

[55] Maitra A 2021 Floating flocks: 2D long-range uniaxial order in 3D active fluids (arXiv:2110.15633)

[56] Ranft J, Basan M, Elgeti J, Joanny J-F, Prost J and Julicher F 2010 Fluidization of tissues by cell division and apoptosis *Proc. Natl Acad. Sci.* **107** 20863

[57] Williams M L and Solnica-Krezel L 2020 Nodal and planar cell polarity signaling cooperate to regulate zebrafish convergence & extension gastrulation movements *eLife* **9** e54445

[58] Merkel M 2014 *PhD Thesis* Technische Universität Dresden

[59] Chien Y H, Keller R, Kintner C and Shook D R 2015 Mechanical strain determines the axis of planar polarity in ciliated epithelia *Curr. Biol.* **25** 2774

[60] Aw W Y, Heck B W, Joyce B and Devenport D 2016 Devenport, transient tissue-scale deformation coordinates alignment of planar cell polarity junctions in the Mammalian skin *Curr. Biol.* **26** 2090

[61] Blanch-Mercader C and Casademunt J 2017 Hydrodynamic instabilities, waves and turbulence in spreading epithelia *Soft Matter* **13** 6913

[62] Bowick M J, Fakhri N, Marchetti M C and Ramaswamy S 2022 Symmetry, thermodynamics, and topology in active matter *Phys. Rev. X* **12** 10501