# Stabilization of active tissue deformation by a dynamic signaling gradient

Muhamet Ibrahimi, Matthias Merkel

► To cite this version:

Muhamet Ibrahimi, Matthias Merkel. Stabilization of active tissue deformation by a dynamic signaling gradient. PRX Life, 2025, 3 (4), pp.043013. 10.1103/yq1s-p2zl. hal-04858096v2

HAL Id: hal-04858096

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

Submitted on 20 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Stabilization of Active Tissue Deformation by a Dynamic Signaling Gradient

Muhamet Ibrahimi<sup>1,2</sup> and Matthias Merkel<sup>1,*</sup>

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

<sup>2</sup>Laboratory of Artificial and Natural Evolution (LANE), Department of Genetics and Evolution, University of Geneva, 1211 Geneva, Switzerland

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

(Received 21 March 2025; accepted 17 October 2025; published 19 November 2025)

A key process during animal morphogenesis is oriented tissue deformation, which is often driven by internally generated active stresses. Yet, such active oriented materials are prone to well-known instabilities, raising the question of how oriented tissue deformation can be robust during morphogenesis. Here we study under which conditions active oriented deformation can be stabilized by the concentration pattern of a signaling molecule, which is secreted by a localized source region, diffuses across the tissue, and degrades. Consistent with earlier results, we find that oriented tissue deformation is always unstable in the gradient-contractile case, i.e., when active stresses act to contract the tissue along the direction of the signaling gradient, and we show that this is true even in the limit of large diffusion. However, active deformation can be stabilized in the gradient-extensile case, i.e., when active stresses act to extend the tissue along the direction of the signaling gradient. Specifically, we show that gradient-extensile systems can be stable when the tissue is already elongated in the direction of the gradient. We moreover point out the existence of a formerly unknown, additional instability of the tissue shape change. This instability results from the interplay of active tissue shear and signal diffusion, and it indicates that some additional feedback mechanism may be required to control the target tissue shape. Taken together, our theoretical results provide quantitative criteria for robust active tissue deformation and explain the lack of gradient-contractile systems in the biological literature, suggesting that the active matter instability acts as an evolutionary selection criterion.

DOI: 10.1103/yq1s-p2z1

## I. INTRODUCTION

Oriented tissue deformation is a key process during animal morphogenesis [1], including, for instance, body axis elongation [2–4] and organ formation [5–10]. In many cases, such tissue deformation is at least in part driven by internally generated active stresses [8,10–13], and in recent years, deforming tissues have successfully been described as active materials [8,14–21]. Moreover, to achieve persistent large-scale oriented deformation, tissues need some kind of tissue-wide orientational information, which can be represented for instance by a polar or a nematic field [8,15,19,22,23]. However, active polar or nematic materials are subject to the well-known Simha-Ramaswamy instability, where an initially globally ordered orientational field loses this order due to active flows [24–26]. This raises the question of how oriented tissue deformation can be stable during animal development.

In animals, oriented tissue deformation can be affected by tissue-scale concentration patterns of signaling molecules. Such patterns of signaling molecules are known to coordinate developmental processes on the tissue or even organism scale [1]. Signaling molecules can provide positional information,

in which case they are called morphogens [3,27–31]. Yet, in some cases, concentration gradients of signaling molecules also provide orientational information that defines the axis of active stresses driving oriented tissue deformation [5–7,12,32], typically by guiding cellular polarity and cell rearrangements [2,4,6,12,33]. Such alignment of orientational information with signaling gradients has recently also been discussed theoretically [34].

In a recent publication [35], we theoretically showed that oriented tissue deformation can be stabilized by signaling gradients. We used a hydrodynamic model for an active polar material, and showed that the gradient of a scalar signaling field can stabilize the otherwise unstable material. We showed that if the scalar field advects with the material flows [36,37], the stability of oriented material deformation depends on the coupling between scalar field gradient direction and deformation axis. In the gradient-extensile case, i.e., when anisotropic active stresses act to extend the material along the gradient direction, the scalar field can stabilize the active material deformation. However, in the gradient-contractile case, i.e., when anisotropic active stresses act to contract the material along the gradient direction, the deformation is unstable. Yet one could always stabilize even a gradient-contractile system by sufficiently increasing the diffusion of the signaling molecule.

In our past publication [35], we have created a uniform scalar field gradient in a simple way through appropriate boundary conditions. While providing advantages for analytical calculations, this is quite different from the real mechanism creating signaling gradients in developing biological tissues [1,28,29,38–40]. Signaling molecules are typically

<footnote>*Contact author: matthias.merkel@cnrs.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Hydrodynamic model. A spatially varying source field $s(r, t)$ secretes a signaling molecule with concentration $c(r, t)$, which also diffuses and degrades. The gradient of $c$ creates extensile/contractile active anisotropic stresses, which drive viscous flows with velocity $v(r, t)$. These flows affect both source and signaling fields through advection. (b) Typical orders of magnitude for tissue deformation rates $\tilde{v}_{xx}^0$ (top), signal diffusion rates $D/L^2$, where $D$ is the diffusion coefficient and $L$ is the linear dimension of the corresponding tissue (middle), and signal degradation rates $k_d$ (bottom) from the biological literature. Blue marks are data from the fruit fly *Drosophila melanogaster*, and orange marks are data from zebrafish, *Danio rerio*. Tissue deformation rates in the fruit fly were estimated from image time series of egg chamber [46,47], dorsal thorax [12], pupal wing [8], hindgut [5,48], renal tubules [7], and germ band [11] and in zebrafish from mesodermal explants [49] and the tailbud [50]. Signaling molecules are Decapentaplegic in the fruit fly [38,40,51] and Cyclops, Squint, Lefty1, and Lefty2 in zebrafish [29,42].</figcaption></figure>

produced by some source region, diffuse across the tissue, and degrade [40,41]. As a consequence, depending on the source region, this can lead to approximately exponentially decaying [38,39,42,43] or periodic [44,45] concentration profiles. However, it is not known how such signaling molecule dynamics would affect stability of active oriented tissue deformation.

In this article, we study the stability of active tissue deformation when accounting for more realistic signaling dynamics [Fig. 1(a)]. We consider an active viscous tissue with local flow velocity $v(r, t)$, where active stresses are defined by the

gradient of a signaling molecule concentration field $c(r, t)$. The signaling field $c(r, t)$ is in turn secreted by a source region that is defined by a field $s(r, t)$, which is conserved and only advects with tissue flows. The signaling field $c(r, t)$ not only advects with tissue flows, but also diffuses and degrades at constant rates everywhere.

Using linear stability analysis, we find that gradient-contractile tissues are always unstable, consistent with our earlier findings [35]. Yet, we also show that for our more realistic signaling dynamics, a gradient-contractile system can never be stabilized by diffusion. We further show that gradient-extensile systems can be stable, in particular if the tissue is longer in the direction of the gradient than transversal to it. Finally, we find that an additional instability can appear for both gradient-extensile and gradient-contractile cases, arising from the interaction between overall tissue shear and diffusion. This instability may point to the necessity of an additional control mechanism for the final tissue dimensions.

In the following, we first introduce the model in Sec. II. We then study the case of a uniform gradient in both $s(r, t)$ and $c(r, t)$ for fixed system size in Sec. III, for which the linear stability can be computed analytically. Next, in Sec. IV, still for fixed system size, we study the more realistic case of a localized source $s(r, t)$, which gives rise to a nonuniform signaling gradient $c(r, t)$. We then discuss the stability of freely deforming systems in Sec. V. Finally, in Sec. VI, we discuss our findings in the context of the biological literature.

## II. MODEL

### A. Hydrodynamic fields and bulk dynamics

We describe the tissue as a 2D active material, where we consider three interacting hydrodynamic fields: a scalar field $s(r, t)$, which describes the source, i.e., cells in the tissue which secrete the signaling molecule; a scalar field $c(r, t)$, which describes the signaling molecule concentration; and a vector field $v(r, t)$, which describes the tissue flows [Fig. 1(a)]. We find that literature values for the respective rates of tissue deformation $\tilde{v}_{xx}^0$, signaling molecule diffusion $D/L^2$, where $D$ is the diffusion coefficient and $L$ denotes the linear dimension of the tissue, and signaling molecule degradation $k_d$ are of roughly similar orders of magnitude [Fig. 1(b)]. This suggests that none of them can be easily neglected, and we thus include all three effects in our model.

The source field $s(r, t)$, which describes signaling molecule secretion, only advects with material flows:

$$ \frac{\partial s}{\partial t} + \partial_i (v_i s) = 0, \quad (1) $$

where $\partial_i$ denotes the partial derivative with respect to the spatial dimension $i \in \{x, y\}$. Here and in the following, we label spatial dimensions by Latin indices $i, j, \dots$ and adopt the Einstein notation of summing over repeated indices.

The signaling field $c(r, t)$ follows a generalized conservation equation:

$$ \frac{\partial c}{\partial t} + \partial_i (v_i c) = k_p s - k_d c + D \partial_i^2 c. \quad (2) $$

Here we have introduced production rate $k_p$, degradation rate $k_d$, and diffusion coefficient $D$. Note that there are several
intra- and extra-cellular mechanisms of signaling molecule transport [29]; here we describe the emergent effective transport on large length scales by advection and diffusion.

Tissue flows $v(r, t)$ are governed by overdamped, incompressible, viscous dynamics:

$$ \partial_i \sigma_{ij} = 0, \quad (3) $$

$$ \partial_i v_i = 0, \quad (4) $$

where the stress tensor is given by

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

Here $\eta$ is the shear viscosity, $\tilde{v}_{ij} := (\partial_i v_j + \partial_j v_i - \partial_k v_k \delta_{ij})/2$ is the shear rate tensor, and $\Pi$ is the hydrostatic pressure, and we use the following expression for the active stress $\tilde{\sigma}_{ij}^a$:

$$ \tilde{\sigma}_{ij}^a = \alpha \left[ (\partial_i c)(\partial_j c) - \frac{1}{2} (\partial_k c)^2 \delta_{ij} \right]. \quad (6) $$

This expression for the active anisotropic stress is the same as in Active Model H [52,53]. Its direction is determined by the sign of the coefficient $\alpha$: for $\alpha < 0$, the activity is gradient-extensile, i.e., it acts to extend the material along the local $c$ gradient direction and contract perpendicularly, while for $\alpha > 0$, the activity is gradient-contractile, i.e., it acts to contract the material along the $c$ gradient direction and extend perpendicularly [Fig. 1(a)].

The expression for the active stress, Eq. (6), is motivated by what is known about developing biological tissues. Specifically, anisotropic internal stresses are usually generated by an anisotropic distribution of cytoskeletal elements within cells. Such an anisotropic distribution is created by cell polarity, which in turn can be controlled by signaling gradients [7,12,33,54,55]. Rather than explicitly including cell polarity as polar field [35], we assume here for simplicity a direct coupling of active stresses to the signaling gradient. This corresponds to the limit where the coupling of active stresses to the signaling gradient is fast compared to the other timescales of the system. While this is clearly a simplification, we do not expect the fundamental results of the analysis to differ much in this limit, as suggested by our earlier work [35]. Note that Eq. (6) represents the lowest-order term for a coupling of an active, anisotropic stress to a scalar field [52].

## B. Boundary conditions

We use periodic boundary conditions, where the periodic box has generally time-dependent dimensions $L_x(t) \times L_y(t)$ and constant area, $L_x(t)L_y(t) = \text{const}$. Throughout this article, we always use one of the following two boundary conditions:

(1) We fix the system dimensions $L_x$ and $L_y$ or

(2) We allow the system to freely deform by fixing the externally applied stress to zero, which results in a pure shear deformation of the periodic box. Integration of Eq. (5) implies that the instantaneous box shear rate is

$$ \frac{1}{L_x} \frac{dL_x}{dt} = - \frac{\alpha}{4\eta L_x L_y} \int \left[ (\partial_x c)^2 - (\partial_y c)^2 \right] dA, \quad (7) $$

where the integral is over the whole periodic box.

We apply fixed system dimensions in Secs. III and IV, and study a freely deforming system in Sec. V.

We investigate the stability of two stationary states: (i) a state where a linear source profile $s$ gives rise to a linear

signaling profile $c$ and (ii) a state where a localized source $s$ generates a spatially decaying signaling field $c$. In both cases, the amplitude of the source field is defined by a constant $s_b$. In case (i), which we study in Sec. III, the boundary conditions are periodic with the exception of the vertical boundaries for the scalar fields $s$ and $c$, where we set for all $y \in [0, L_y)$:

$$ s(0, y) = s(L_x, y) - s_b, \quad (8) $$

$$ c(0, y) = c(L_x, y) - \frac{k_p}{k_d} s_b, \quad (9) $$

with fixed $s_b$. These modified conditions remove the discontinuity that would otherwise appear for linear profiles in $s$ and $c$ at the boundary, and thus ensure that linear $s$ and $c$ profiles can be stationary [35]. In case (ii), which we study in Secs. IV and V, we use standard periodic boundary conditions without offset.

## C. Dimensionless equations

We nondimensionalize the set of equations Eqs. (1)–(9) in the following way:

(1) We choose the unit of length such that the system area is one, $L_x L_y = 1$.

(2) We choose the unit of the signaling field $c$ such that $k_p = k_d$.

(3) We choose the units of time and of the source field $s$, such that the free deformation rate of a squared system is 1. To this end, we set $|\alpha|/4\eta = 1$, and in case (i) of a linear gradient, we also set $s_b = 1$, while in case (ii) of a nonuniform gradient, we fix $s_b$ as described in Sec. IV (Appendix E 1).

The last condition allows to compare the linear stability of systems with different source profiles, because the reference timescale is always the respective free deformation rate.

Rescaled accordingly, the dimensionless equations read

$$ \frac{\partial s}{\partial t} + v_i \partial_i s = 0, \quad (10) $$

$$ \frac{\partial c}{\partial t} + v_i \partial_i c = (D \partial_i^2 - k_d) c + k_d s, \quad (11) $$

$$ 0 = \frac{1}{4} \partial_i^2 v_j - \partial_j \Pi' + \operatorname{sgn}(\alpha) \partial_i [(\partial_i c)(\partial_j c)], \quad (12) $$

$$ \partial_i v_i = 0. \quad (13) $$

Here we have introduced the sign function, $\operatorname{sgn}(\alpha) := \alpha/|\alpha|$, and have set $\Pi' := \Pi + \operatorname{sgn}(\alpha)(\partial_k c)^2/2$. The offset introduced at the boundary conditions in Eqs. (8) and (9) is rescaled to 1 for both scalar fields.

## III. LINEAR SOURCE PROFILE

To build intuition, we study for fixed box dimensions $L_x = L_y = 1$ the linear stability of a stationary state with linear spatial profiles in $s$ and $c$ [Fig. 2(a)], given by

$$ s_0 = x, \quad c_0 = x, \quad \mathbf{v}_0 = \mathbf{0}. \quad (14) $$

There are no flows in the stationary state, because the gradient of $c$ is homogeneous, and thus the active stress is homogeneous across the system of fixed dimensions.
<figure><img src="image_4.png" /><figcaption>FIG. 2. Stability of a fixed-size system with a linear source profile: gradient-extensile systems are marginally stable, whereas gradient-contractile systems are unstable. (a) The linear source profile, $s_0(x) = x$ in dimensionless units. In the stationary state, $c_0(x) = s_0(x)$. (b) Maximal perturbation growth rates $\omega_{\max}$ in the $D-k_d$ parameter space show that gradient-extensile systems are marginally stable, i.e., $\omega_{\max} = 0$, highlighted in light gray. (c) $\omega_{\max}$ as a function of signaling molecule diffusion $D$ for a gradient-contractile system with $k_d = 1$, plotted in log-log scaling. (d) $\omega_{\max}$ in the $D-k_d$ parameter space show that gradient-contractile systems are always unstable. The dashed gray line corresponds to the $k_d = 1$ curve plotted in panel (c). Panels (b)–(d) were created from the result in Appendix D2, where the perturbation wave vector $q_{\min} = (0, 2\pi)$ in dimensionless units was used to maximize the perturbation growth rate.</figcaption></figure>

In this case, linear stability can be computed analytically. When we add a small perturbation around this state,

$$s = x + \delta s, \quad c = x + \delta c, \quad v = \delta v, \qquad (15)$$

and linearize the dynamics, only constant prefactors occur. As a consequence the solutions in $(\delta s, \delta c, \delta v)$ are spatial Fourier modes with wave vectors $q = q(\cos \phi, \sin \phi)$, whose amplitudes grow or shrink exponentially with time $t$ (Appendix D).

We first discuss the limit without diffusion, $D = 0$. When we subtract the linearized versions of source and signaling dynamics, Eqs. (10) and (11), we obtain (Appendix C)

$$\frac{\partial}{\partial t}(\delta c - \delta s) = -k_d(\delta c - \delta s). \qquad (16)$$

Thus, any difference between signaling field $c$ and source field $s$ decays at rate $k_d$. In addition, the combination of signaling dynamics and generated active flows, Eqs. (11)–(13), leads to the same results as for a single scalar field discussed in earlier publications [35,53] (Appendix D1). Specifically, for the combined perturbation field $\delta c_{\text{eff}} := \delta c + k_d/\omega_{\text{act}}\delta s$, we obtain [Eq. (D8) in Appendix D]

$$\frac{\partial \delta c_{\text{eff}}}{\partial t} = \omega_{\text{act}} \delta c_{\text{eff}}, \qquad (17)$$

where $\omega_{\text{act}} := 4 \operatorname{sgn}(\alpha) \sin^2 \phi$ is the perturbation growth rate created by the active flows. It takes the same form with-

out source field [35]. Thus, without diffusion, the system is marginally stable in the gradient-extensile case, i.e., the maximal perturbation growth rate is zero, and unstable in the gradient-contractile case.

In the presence of diffusion, $D > 0$, gradient-extensile systems remain marginally stable [Fig. 2(b)]. However, gradient-contractile systems can not be stabilized by diffusion, even though the maximal perturbation growth rate decreases with diffusion [Figs. 2(c) and 2(d)]. This is in contrast to earlier findings for a system with a boundary-imposed signaling gradient, where we found that diffusion can stabilize gradient-contractile systems [35]. To understand where the difference comes from, we consider the limit of very large $D$ (the general case is discussed in Appendix D). The linearized version of the signaling dynamics, Eq. (11), in Fourier space reads

$$\frac{\partial \delta c}{\partial t} + \delta v_x = -(Dq^2 + k_d)\delta c + k_d \delta s. \qquad (18)$$

The flow field $\delta v$ enters on the left-hand side due to advection, where $\partial_x c_0 = 1$ in dimensionless units. The flow $\delta v$ is actively created by the signaling field itself through Eqs. (12) and (13), and its $x$ component in Fourier space reads $\delta v_x = -\omega_{\text{act}}\delta c$, with the active flow timescale $\omega_{\text{act}}$ introduced before. Denoting the minimal nonzero wave vector magnitude in our system by $q_{\min}$, we take the limit where $Dq_{\min}^2$ is much larger than all other relevant timescales in the system. In this case, according to Eq. (18), the Fourier amplitude $\delta c$ for each wave vector relaxes adiabatically fast towards

$$\delta c \approx \frac{k_d}{Dq^2} \delta s. \qquad (19)$$

This effectively corresponds to a system where we again just have a single scalar field $s$, which advects with active tissue flows $v$, but does not diffuse. The active tissue flows are generated by the perturbations of $c$, and their $x$ component is given by $\delta v_x = -\omega_{\text{act}}\delta c \approx -(\omega_{\text{act}}k_d/Dq^2)\delta s$. Inserting this into the linearized version of the source dynamics in Fourier space, we obtain [compare to Eq. (17)]

$$\frac{\partial \delta s}{\partial t} = -\delta v_x = \frac{\omega_{\text{act}} k_d}{Dq^2} \delta s. \qquad (20)$$

Thus, the perturbation growth rate is given by the growth rate without diffusion, $\omega_{\text{act}}$, multiplied by $k_d/Dq^2$ [Fig. 2(c)]. As a consequence, gradient-contractile systems are unstable even in the limit of a large diffusion constant. The reason for this is that the source field itself advects with the active flows but does not diffuse. For the same reason, gradient-extensile systems remain only marginally stable even in the presence of signaling molecule diffusion.

## IV. LOCALIZED SOURCE PROFILE

We next discuss how the linear stability changes when the source profile is not linear but is instead localized to some region of finite width $w$. This is motivated by what we know about developing biological tissues [1,38,41,45]. In this section, we discuss a fixed-size system $L_x = L_y = 1$, while in the next one (Sec. V) we discuss a freely deforming system.

As the initial state, we use a source profile defined by a von Mises distribution, which is a generalization of a Gaussian to
<figure><img src="image_5.png" /><figcaption>FIG. 3. Stability of a fixed-size system with a localized source profile. (a) The localized source profile $s_0(x)$ from Eq. (21) for different values of the source width $w$. (b) The localized signaling profile for different values of diffusion $D$. Here $w = 1/2\pi$ in all cases. (c), (d) Snapshots of the source profile from numerical simulations of the perturbed stationary state (with $w = 1/2\pi$, $D = 10^{-2}$, $k_d = 1$). While the source profile appears to be maintained in the (c) gradient-extensile case, it gets distorted in the (d) gradient-contractile case. (e) Effect of aspect ratio $L_x/L_y$ on maximal perturbation growth rates $\omega_{max}$ for different source widths. Gradient-contractile systems (orange) are always unstable, while gradient-extensile systems (green) have in some regions an $\omega_{max}$ that can numerically not be distinguished from zero. Parameter values: $k_d = 1$ and $D = 1$. (f) Effect of the source width $w/L_x$ on maximal perturbation growth rate for different box aspect ratios. Parameter values: $k_d = 1$ and $D = 1$. (g, h) Color plots of maximal perturbation growth rates for gradient-extensile (g) and gradient-contractile (h) couplings with varying source widths (horizontal axis) and aspect ratios (vertical axis). Marginally stable regions ($\omega_{max} = 0$) are highlighted in light gray. Solid black curves represent contours, whereas dashed gray lines corresponds to solid curves in panels (e) and (f). Parameter values: $k_d = 1$ and $D = 1$. (i), (j) Maximal perturbation growth rates in the $D-k_d$ parameter space for gradient-extensile (i) and gradient-contractile (j) couplings. Numerically marginally stable regions are highlighted in light gray. Parameter values: $L_x/L_y = 1$ and $w/L_x = 0.05$.</figcaption></figure>

periodic boundary conditions:

$$s_0(\mathbf{r}) = s_b \left[ \exp\left(\frac{\cos(q_0 x)}{(q_0 w)^2}\right) - \exp\left(-\frac{1}{(q_0 w)^2}\right) \right] \quad (21)$$

with $q_0 := 2\pi/L_x$. Here $w$ is a parameter that adjusts the width of the source region, the second term in Eq. (21) ensures that $s_0(\mathbf{r})$ is zero at its minimum, and $s_b$ is a prefactor whose value is fixed such that the box deformation rate is normalized to one [Eq. (E1) in Appendix E 1]. As initial signaling field,
<figure><img src="image_6.png" /><figcaption>FIG. 4. Instability in a gradient-extensile system with fixed size and a localized source profile. (a) Snapshots of the late-time dynamics of the system in Fig. 3(c), showing the source profile at $t = 150$ (left) and $t = 1500$ (right). (b) Perturbation mode with maximum growth rate for system from Fig. 3(c) and panel (a). The color plot shows the signaling molecule perturbation, $\delta c$, and white arrows show the flow field perturbation, $\delta v$.</figcaption></figure>

we choose the $c$ field that is stationary in the absence of flows for the source given by Eq. (21). In Fourier space,

$$c_0(\boldsymbol{q}) = \frac{k_d}{k_d + Dq^2} s_0(\boldsymbol{q}). \qquad (22)$$

Thus, for $D = 0$, signaling and source profiles coincide, and examples for $s_0(x) = c_0(x)$ for different $w$ values are shown in Fig. 3(a). Meanwhile, finite diffusion $D > 0$ smoothens the signaling profile $c_0(x)$ as compared to the source profile $s_0(x)$ [Fig. 3(b)].

In this initial state, there are no flows in the system, $v_0 = 0$. This is because the active stress in this case is parallel to the $x$ axis and homogeneous along the $y$ axis, because of incompressibility, and because of the fixed system size, which together imply that there are no Stokes flows (Appendix B 2).

Numerical simulations of the full dynamics, Eqs. (10)–(13), are shown in Figs. 3(c) and 3(d) and in Movies S1 and S2 (shown is the source field) [56]. We find that with chosen parameters ($w = 1/2\pi$, $k_d = 1$, and $D = 10^{-2}$), the gradient-extensile system appears to be stable until simulation time $t = 3$ [Fig. 3(c)], yet an instability becomes apparent at much later times [Fig. 4(a) and Movie S1 [56]]. Meanwhile, the gradient-contractile system shows an instability already early on [Fig. 3(d) and Movie S2 [56]].

For a more systematic characterization of the system's stability, we perform again a linear stability analysis. While the eigenmodes of the linearized dynamics still grow or shrink exponentially with time, they can only be represented by a superposition of several spatial Fourier modes (Appendix C 2). We thus compute the perturbation eigenmodes and their respective growth rates numerically (Appendix C 3). Note that due to the numerical Fourier-space convolutions, we also find an aliasing effect arising from a finite spatial discretization, which we treat using standard methods in the field (Appendix C 3).

We find again that gradient-contractile systems are always linearly unstable [Figs. 3(e), 3(f), 3(h), and 3(j)]. Specifically, systems are more unstable when the source width is below $w \lesssim 0.06L_x$ [Figs. 3(f) and 3(h)]. Moreover, increasing the diffusion coefficient $D$ decreases the perturbation growth rate, but it never reaches zero. Instead, in contrast to the linear gradient (Sec. III), the maximum growth rate converges to a finite constant in the limit of large diffusion, $D \to \infty$ (Appendix E 2). It converges to a finite constant, because we nondimensionalize rates with respect to the free deformation rate, which also decreases as diffusion $D$ becomes larger due to an increased smoothing of $c$ (Appendix E 2). This is different from the linear-gradient scenario (Sec. III), where a constant signaling gradient and thus a constant free deformation rate are imposed by the boundary conditions.

For gradient-extensile systems, the situation is more complex than for a linear gradient. Specifically, we find that the stability depends on the aspect ratio of the system. For aspect ratios $L_x/L_y \lesssim 2...3$, the system can be unstable, while it is numerically marginally stable for aspect ratios above [Figs. 3(e) and 3(g)]. Here, by “numerically marginally stable,” we mean that we numerically cannot distinguish between marginally stable and numerically small, but finite growth rates—indeed, this distinction is also irrelevant in practice for a real biological system. For aspect ratios between $\approx 1$ and $2...3$, the system is effectively unstable at intermediate source widths $w/L_x$ between 0.04 and 0.2 and otherwise numerically marginally stable [light gray regions in Fig. 3(g)]. For even lower aspect ratios, we only observe numerical marginal stability for very narrow sources with $w/L_x \lesssim 0.04$ [Figs. 3(f) and 3(g)]. The type of perturbation creating the most unstable mode appears to be qualitatively always the same: recirculating active flows create spatial variations of the source region width, which create corresponding modulations in the signaling field, leading to the recirculating active flows [Figs. 4(a) and 4(b)]. Finally, for given aspect ratio and source width, the system is generally more stable for small signaling length scales $\sqrt{D}/k_d$ [Fig. 3(i)]. In the large-diffusion limit, the maximum perturbation growth rate converges to a constant for the same reason as in the gradient-contractile case (Appendix E 2).

Note that even if gradient-extensile systems are formally unstable in a substantial part of the parameter space, for aspect ratios close to one the maximal perturbation growth rates can be several orders of magnitude smaller as compared to the gradient-contractile case. For instance, the maximal growth rate for the example simulation in Figs. 3(c), 4(a), and Movie S1 [56] is positive, but very small, $\omega_{max} = 1.9 \times 10^{-3}$ [as
compared to $\omega_{max} = 4.2$ in the gradient-contractile example in Fig. 3(d)].

Taken together, the gradient-extensile system is effectively stable when it is more elongated in the direction of the gradient, or when the signaling length scale $\sqrt{D/k_d}$ is small.

## V. FREELY DEFORMING SYSTEM

Here we study the linear stability of a system with a localized source $s_0$, Eq. (21), under free deformation of the periodic box.

Even for homogeneous shear deformation, the source profile $s_0$ is not stationary, but becomes distorted due to advection with the box deformation [Fig. 5(a)]. To address this, following our earlier work [35], we define the codeforming coordinates, $\bar{r} = (\bar{x}, \bar{y})$, based on the laboratory-frame coordinates $r = (x, y)$ by rescaling both axes:

$$ \bar{x} := \frac{L_x(0)}{L_x(t)}x = l_x^{-1}(t)x, \quad (23) $$

$$ \bar{y} := \frac{L_y(0)}{L_y(t)}y = l_y(t)y, \quad (24) $$

where $l_x(t) := L_x(t)/L_x(0)$ is the box shear. Moreover, we also introduce a flow field in codeforming coordinates, $\bar{v}(\bar{r})$, which corresponds to nonaffine flows that may occur in addition to the pure shear deformation of the periodic box (Appendix A 3). Thus, for $\bar{v} = 0$ the system deforms affinely, and the source field $s_0(\bar{r})$ is stationary with respect to the codeforming coordinates [Fig. 5(b)].

Yet, in the general case of a finite diffusion constant, the signaling field $c$ does not reach a stationary state—not even in codeforming coordinates. This is because diffusion as observed in codeforming coordinates depends on box shear $l_x$ through the transformation rules in Eqs. (23) and (24). Thus, as a consequence of the permanently changing $l_x$, there is no stationary state in the freely deforming case. Hence, no standard linear stability analysis is possible—we instead analyze the local Lyapunov exponents of the system.

To compute the local Lyapunov exponents, we consider the system state by combining source and signaling profiles in codeforming coordinates with box shear, $\bar{z}(t) = [s(\bar{r}, t), c(\bar{r}, t), l_x(t)]$. This state constantly evolves through Eqs. (7) and (10)–(13), which we summarize here as $d\bar{z}/dt = \bar{F}(\bar{z})$. As initial condition, we use the source profile $s_0$ from Eq. (21), and the corresponding signaling profile $c_0$ that would be stationary in a fixed system, Eq. (22). To simplify our discussion here, we will focus on the limit $w \to \infty$, where the source $s_0(r)$ given by Eq. (21) corresponds to a cosine profile, $s_0(r) \sim \cos(2\pi x/L_x) + 1$. The corresponding signaling profile is then also a cosine profile. We discuss the unperturbed box dimension dynamics, $l_x(t)$, in Appendix F1, and focus here only on the perturbation dynamics.

Local Lyapunov exponents characterize whether and how quickly two slightly different solutions, $\bar{z}$ and $\bar{z}'$, converge or diverge. The difference between both, $\delta\bar{z} := \bar{z}' - \bar{z}$, follows the dynamics $d\delta\bar{z}/dt = J(\bar{z}) \cdot \delta\bar{z}$, where $J(\bar{z})$ is the Jacobian of $\bar{F}(\bar{z})$. The local Lyapunov exponents at $\bar{z}$ are then given by the eigenvalue spectrum of $J(\bar{z})$. When applied to a stationary state, this approach coincides with linear stability analysis.

<figure><img src="image_7.png" /><figcaption>FIG. 5. Stability of a freely deforming system with a diffusive signaling molecule, studied in codeforming coordinates. (a) In laboratory-frame coordinates, a localized source $s(r, t)$ is not stationary when the system deforms affinely, since its width changes proportionally with the system dimension $L_x$ due to advection. (b) However, in codeforming coordinates, Eqs. (23) and (24), the localized source $s(\bar{r}, t)$ remains stationary under affine deformation. (c) Illustration of the affine diffusion-shear instability: diffusion reduces the amplitude $\hat{c}$ of the signaling field, but a larger $L_x$ decreases the effective diffusion coefficient observed in codeforming coordinates. Thus, a larger $L_x$ leads to a larger $\hat{c}$, which in turn increases active shear and thus $L_x$. (d) Local Lyapunov exponents (LLE) for the gradient-extensile system over time. Blue crosses indicate the numerically computed LLE of the affine diffusion-shear perturbation with the largest growth rate, and the light blue curve displays the analytical prediction according to Eq. (F9) in Appendix F2a. Green squares denote the maximal LLE of all other perturbations. (e) The same results as in panel (d), but for the gradient-contractile system. Specifically, black crosses indicate the LLE of the affine diffusion-shear perturbation with the largest growth rate, and orange squares denote the maximal LLE of all other perturbations.</figcaption></figure>

In the limit of no diffusion, $D = 0$ (previous section, Sec. IV), the linearized dynamics of source and signaling, Eqs. (C4) and (C5) (Appendix C1), do not depend on the
perturbation of the box shear, $\delta l_x$. Thus, the perturbation dynamics of s and c decouples from that of $l_x$, and most of the Lyapunov spectrum corresponds to the eigenvalue spectrum of the fixed system (Sec. IV). The remaining eigenvalue of the local Lyapunov spectrum is given by the system size perturbation dynamics:

$$ \frac{d\delta l_x}{dt} = \frac{\operatorname{sgn}(\alpha)}{l_x^2(t)L_x^2(0)} \delta l_x. \qquad (25) $$

The prefactor in this equation is positive for gradient-contractile systems and negative for gradient-extensile systems. Hence, for $D = 0$, the gradient-contractile system has positive local Lyapunov exponents, i.e., the system dynamics strongly depends on the initial conditions. Meanwhile, for gradient-extensile systems, the maximal local Lyapunov exponent is given by the maximal eigenvalue for the fixed-size system (Sec. IV).

### A. Affine diffusion-shear instability

For finite diffusion, $D > 0$, the equations defining the linearized dynamics, Eqs. (C4)–(C6), are all coupled to each other and the full Jacobian needs to be diagonalized in order to obtain the local Lyapunov exponents. At each time point of the trajectory $\bar{z}(t)$, we plot the maximal Lyapunov exponent in Figs. 5(d) and 5(e).

We first discuss the gradient-extensile case, where we find that initially, perturbation growth rates can be positive [blue crosses and green squares in Fig. 5(d)]. However, at later times, as the system shears and $l_x$ grows sufficiently, all perturbation growth rates become maximally zero [Fig. 5(d)]. Initially, i.e., at low $l_x$, we found that the fastest-growing perturbation mode corresponds to a new kind of instability that we call “affine diffusion-shear instability” [blue crosses in Fig. 5(d)]. The nature of this instability is quite different from the active scalar matter instability [35,53]—specifically, it affects neither the nonaffine part of the flow field $\bar{v}$ nor the source field s. Moreover, the signaling profile keeps its cosine shape, and only its amplitude $\hat{c}$ and the box shear $l_x$ are affected by the perturbation. The instability arises because diffusion diminishes the amplitude of the signaling cosine profile, $\hat{c}$. Thus, increasing $l_x$ by some positive amount $\delta l_x$ decreases the effect of diffusion, which acts to increase $\hat{c}$, which in turn increases the active stresses making $l_x$ grow further [Fig. 5(c)]. In Appendix F2a, Eq. (F9), we obtain an analytical prediction for the growth rate of this instability [light blue curve in Fig. 5(d)], which matches our numerical results (blue crosses). From our analytical prediction, we find that this instability appears as long as the diffusion length $\lambda := \sqrt{D/k_d}$ is larger than $L_x$, or more precisely $\lambda > L_x/(2\pi\sqrt{3})$ (Appendix F2a). Remarkably, this instability does not disappear even when introducing an arbitrary nonlinearity into the active stress (Appendix G).

For the gradient-contractile case, we find several modes with a positive growth rate [orange squares in Fig. 5(e)]. Moreover, we again have an eigenmode of the Jacobian that is an affine diffusion-shear perturbation (black crosses in Fig. 5(e); compare to the predicted growth rate marked by the light blue curve).

Finally, in Appendix F2b we study the affine diffusion-shear instability for the general case of a finite source width w. We find a positive growth rate whenever the order of magnitude of the diffusion length becomes larger than the source width, $\lambda \gtrsim w$ (Appendix F2b).

### VI. DISCUSSION

In developing animals, biological tissues can undergo active anisotropic deformation [1]. However, classical active matter results predict that such deformations are unstable under certain conditions [24,25,53]. In several systems, the local direction of active tissue deformation is controlled by the gradient of a signaling molecule concentration field [7,12,33]. Such signaling molecules are typically produced by the cells in a localized source region of the tissue, and they can form gradients through diffusion and degrade [1]. However, both signaling pattern and source region are also advected by the flows generated by the active deformation [36]. Here we study the stability of a system which includes source field, signaling field, and signaling-gradient-controlled active tissue deformation. We find that systems where the active stresses are gradient-contractile are always unstable, while gradient-extensile systems can be marginally stable under certain conditions. Our findings are consistent with, and potentially explain, an observed abundance of gradient-extensile systems in the biological literature [3,5,7,32,44,47,49,54,57,58] and the virtual absence of gradient-contractile systems (also discussed in [59]), suggesting that the instability discussed here may act as an evolutionary selection criterion.

In our previous work, we studied a simplified signaling dynamics, where the overall gradient was imposed by the boundary conditions [35]. Here, we show that more realistic signaling dynamics, including secretion, diffusion, and degradation, lead to a number of important differences. First, in the gradient-contractile case, neither degradation nor diffusion can stabilize the system. While diffusion can decrease the growth rates of perturbations, we find that it never suffices to fully stabilize the system. This is because both signal and source fields are advected with the active flows, while only the signal diffuses. In the case of a linear source profile, the maximal perturbation growth rate could become arbitrarily small for a sufficiently large diffusion coefficient. However, in the case of a localized source, the maximal perturbation growth rate converges to a finite multiple of the active shear rate in the limit of a large diffusion coefficient (Appendix E2).

Second, we show that for a localized source region, gradient-extensile systems can become unstable whenever the tissue is too short in the direction of the gradient as compared to the transversal direction. Conversely, we numerically obtain marginal stability whenever the tissue is more elongated in the direction of the gradient or when the signaling length scale, $\sqrt{D/k_d}$, is small. However, we also find that even in the unstable regime, the maximum growth rate of gradient-extensile systems are often orders of magnitude smaller than the active shear rate, in particular for aspect ratios close to one [Figs. 3(e)–3(j)]. It is not impossible that during development, biological tissues transiently accept an instability with a small growth rate until entering a stable regime once they become more elongated.
Finally, in deforming systems, an “affine diffusion-shear instability” can appear. Different from the other instabilities discussed here, which are all related to the Simha-Ramaswamy instability, the affine diffusion-shear instability does not require nonaffine tissue flows. Instead, it is only reflected in an interplay between the amplitude of the signaling molecule and the overall tissue shear rate. This instability appears in tissues with a localized source whose width is smaller than the signaling diffusion length $\sqrt{D/k_d}$. It occurs even when adding any nonlinearity in the dependence of active stress on the signaling gradient (Appendix G). In many developing tissues, the signaling profile has a larger extent than the source region [1,5,38], and our work suggests that this instability should occur. Although this instability does not perturb the patterns of source, signaling molecule, or flow, it perturbs their amplitude and, most importantly, leads to a lack of control of the final tissue dimensions. This points to an additional feedback mechanism not included in our model that may be required to control the final tissue shape.

Our study further motivates experiments, and a key question is whether there are any gradient-contractile biological systems at all. There are many systems where signaling molecules are known to be necessary for active anisotropic tissue deformation. However, in most systems, it is not yet known whether the signal controls the directionality or only the magnitude of the active stresses. For instance, there are still open debates about what controls the directionality of active stresses in what is probably the best-studied anisotropic tissue deformation process, which is germ band extension (GBE) in the fruit fly *Drosophila melanogaster*. Recent work on GBE proposed that a dorso-ventral signaling gradient modulates the magnitude of junctional myosin, which plays the role of active stresses [60]. Yet, different mechanisms have been proposed for the control of the myosin anisotropy and its directionality. Classically, these include anterior-posterior-oriented gradients [13,44], which would correspond to a gradient-extensile coupling. But more recently, a group has proposed that a dorso-ventrally aligned tissue tension could provide a global mechanical orienting signal [36,60]. Thus, in many cases, more experiments are needed to cleanly disentangle the signals that merely modulate the amplitude of active stresses from those that control their direction. For instance, to modify the direction of the signaling pattern while keeping a similar overall magnitude, one could combine null-mutants with a graded over-expression [33] or the external supply of a gradient, e.g. through signaling-molecule-soaked beads [1,27]. We believe that systematically combining insights from active matter theory with experiments in this way will help decode the biological mechanisms underlying oriented tissue deformation and its robustness.

## ACKNOWLEDGMENTS

We thank the Centre Interdisciplinaire de Nanoscience de Marseille (CINaM) for providing office space. 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. This project was also supported by

the grant RobustTissue attributed to M.M. by the French National Research Agency (ANR-22-CE30-0039).

## DATA AVAILABILITY

The code that supports the findings of this article is openly available [61].

# APPENDIX A: CODEFORMING COORDINATES

## 1. Definition

Following Ref. [35], we introduce the codeforming coordinate system, $(\bar{\mathbf{r}}, \bar{\mathbf{t}}) = (\bar{x}, \bar{y}, \bar{\mathbf{t}})$, to analytically solve the linearized dynamics of systems under pure shear deformation.

The codeforming coordinates map to the laboratory coordinates $(\mathbf{r}, t) = (x, y, t)$ in the following way:

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

$$t = \bar{\mathbf{t}}, \quad (A2)$$

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

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

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

Note that for fixed system dimensions, where $L_x(t) = L_x(0)$, and thus $l_x(t) = 1$, both, codeforming and laboratory coordinates coincide, $\bar{\mathbf{r}} \equiv \mathbf{r}$.

## 2. Partial derivatives

As a direct consequence of Eqs. (A1) and (A2), partial derivatives of some quantity $f$ transform as

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

$$\bar{\partial}_t f := \frac{\partial f(\bar{\mathbf{r}}, \bar{\mathbf{t}})}{\partial \bar{\mathbf{t}}} = \partial_t f + (\partial_i f) \dot{s}_{ij} \bar{r}_j, \quad (A5)$$

where $\partial_i f := \partial f(\mathbf{r}, t)/\partial r_i$, $\partial_t f := \partial f(\mathbf{r}, t)/\partial t$, and $\dot{s}_{ij} := ds_{ij}/dt = ds_{ij}/d\bar{\mathbf{t}}$. Thus, the partial time derivative in codeforming coordinates, $\partial_t f$, 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 laboratory coordinates.

## 3. Velocity and velocity gradient

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 Eq. (A1):

$$v_i = \dot{s}_{ij} \bar{r}_j + s_{ij} \bar{v}_j, \quad (A6)$$

where $\bar{v}_i := d\bar{r}_i/dt = d\bar{r}_i/d\bar{\mathbf{t}}$ is the codeforming velocity, with $\bar{r}(\bar{\mathbf{t}})$ being the tracer trajectory in codeforming coordinates.

The first term in Eq. (A6) corresponds to a motion due to the affine transformation according to box coordinates. Thus, $\bar{v}_i$ can be interpreted as the nonaffine component of the flow
field. This can also be seen more explicitly by computing the velocity gradient from Eq. (A6):

$$ \partial_i v_j = v_{ij}^0 + s_{li}^{-1} s_{jk} \bar{\partial}_l \bar{v}_k, \quad (A7) $$

where $v_{ij}^0 := s_{ki}^{-1} s_{jk}$ is the average box deformation rate tensor. For $s_{ij}$ as defined in Eq. (A3), the box deformation rate tensor is

$$ v_{ij}^0 = \begin{pmatrix} \dot{l}_x/l_x & 0 \\ 0 & -\dot{l}_x/l_x \end{pmatrix} = \begin{pmatrix} \dot{L}_x/L_x & 0 \\ 0 & -\dot{L}_x/L_x \end{pmatrix}. \quad (A8) $$

### 4. Total derivative

To obtain a transformation formula for the convective derivative, we consider again our tracer and the presence of some spatiotemporal field $f$. The convective derivative corresponds to the total derivative of the value of $f$ that the tracer locally sees. Thus, we expect analogous expressions for the convective derivative in both laboratory and codeforming systems, $\dot{f} := df/dt = df/d\bar{r}$. Indeed, using Eqs. (A4)–(A6), we obtain

$$ \dot{f} = \partial_t f + v_i(\partial_i f) = \bar{\partial}_t f + \bar{v}_i(\bar{\partial}_i f). \quad (A9) $$

### 5. Fourier transform

Moreover, we define the codeforming Fourier transformation of a quantity $f$ such that

$$ f(\bar{r}) = \sum_{\bar{q}} f(\bar{q}) e^{i\bar{q}\cdot\bar{r}}, \quad (A10) $$

where the sum is over all codeforming wave vectors $\bar{q} = (\bar{q}_x, \bar{q}_y)$ with $\bar{q}_i = 2\pi n_i/L_i(0)$, where $n_i \in \mathbb{Z}$ and $i \in \{x, y\}$. A codeforming Fourier mode with wave vector $\bar{q}$ is distorted over time by the overall system deformation:

$$ q_x(\bar{q}, t) = l_x^{-1}(t)\bar{q}_x, \quad (A11) $$

$$ q_y(\bar{q}, t) = l_x(t)\bar{q}_y. \quad (A12) $$

As in Eq. (A10), we have the usual derivation rule, where the Fourier transform of $\bar{\partial}_j f(\bar{r}, \bar{r})$ is $i\bar{k}_j f(\bar{q}, \bar{r})$. From Eqs. (A10) and (A1), it also follows that a given codeforming Fourier mode with wave vector $\bar{q}$ corresponds to a laboratory-frame Fourier mode with wave vector $\bar{q}$ with components

$$ q_i = \bar{q}_j s_{ji}^{-1}, \quad (A13) $$

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

## APPENDIX B: DYNAMICS IN CODEFORMING COORDINATES

Using Appendix A, we can express the dimensionless dynamic equations, Eqs. (10)–(13), in codeforming coordinates:

$$ \bar{\partial}_t s + \bar{v}_i(\bar{\partial}_i s) = 0, \quad (B1) $$

$$ \bar{\partial}_t c + \bar{v}_i(\bar{\partial}_i c) = (D_s \bar{\partial}_i^2 \bar{\partial}_j - k_d) c + k_d s, \quad (B2) $$

$$ 0 = \frac{1}{4} s_{ij} \partial_i^2 \bar{v}_l - \partial_j \Pi + \partial_i \bar{\sigma}_{ij}^a, \quad (B3) $$

$$ s_{il} \partial_i \bar{v}_l = 0, \quad (B4) $$

where $\bar{\sigma}_{ij}^a$ is given by Eq. (6). Equations (B3) and (B4) are still expressed using laboratory-frame derivatives, because this will allow for a more convenient solution of the nonaffine flows in the next section.

Finally, the system size dynamics for the freely deforming system, Eq. (7), becomes

$$ \frac{dl_x}{dt} = -\operatorname{sgn}(\alpha) l_x \iint [l_x^{-2}(\bar{\partial}_x c)^2 - l_x^2(\bar{\partial}_y c)^2] d\bar{x} d\bar{y}, \quad (B5) $$

where the integral is over the box at time point 0, i.e., the integration variables are the codeforming coordinates $\bar{x}$ and $\bar{y}$.

### 1. Nonaffine flows

For deforming systems, we will generally consider the dynamics of source and signaling fields in codeforming coordinates, Eqs. (B1) and (B2). To obtain the nonaffine flows, we solve Eqs. (B3) and (B4) by taking the laboratory-frame Fourier transform, which yields after some transformations

$$ \bar{v}_k = \frac{4i s_{ki}^{-1} q_l}{q^2} \left( \delta_{ij} - \frac{q_l q_j}{q^2} \right) \bar{\sigma}_{lj}^a \quad (B6) $$

with $\bar{\sigma}_{ij}^a$ given by Eq. (6).

Note that Eq. (B6) is undefined for $q = 0$. Indeed, any constant homogeneous velocity can always be added to the flow field. Here and in the following, we will ignore this, because it does not affect the physics in any way.

### 2. Absence of nonaffine flows for $\partial_y c = 0$

If the signaling field $c$ is homogeneous in $y$ direction, i.e., $\partial_y c = \partial_{\bar{y}} c = 0$, then there are no nonaffine flows, i.e., there are no flows in codeforming coordinates, $\bar{v} = 0$. For fixed system dimensions, this corresponds to a complete absence of flows, $\bar{v} = 0$.

Formally, the absence of nonaffine flows follows from Eq. (B6), which links $\bar{v}$ in Fourier space to the stress tensor in Fourier space. From Eq. (6) we obtain for $\partial_y c = 0$, that $\bar{\sigma}_{xy}^a = 0$ everywhere, and $\bar{\sigma}_{xx}^a$ is independent of $y$. Thus, in Fourier space, $\bar{\sigma}_{xx}^a$ is nonzero only for wave vectors $\bar{q}$ that have a vanishing $y$ component, $q_y = 0$. Insertion in Eq. (B6) yields indeed $\bar{v} = 0$.

## APPENDIX C: LINEARIZED DYNAMICS

### 1. In codeforming real space coordinates

We perturb the system around the dynamic state $\bar{z}_0 = (s_0(\bar{r}), c_0(\bar{r}), l_{x0})$ to linear order, where $s_0$ and $c_0$ only depend on $\bar{x}$:

$$ s(\bar{r}) = s_0(\bar{x}) + \delta s(\bar{r}), \quad (C1) $$

$$ c(\bar{r}) = c_0(\bar{x}) + \delta c(\bar{r}), \quad (C2) $$

$$ l_x = l_{x0} + \delta l_x. \quad (C3) $$

Because $c_0$ only depends on $\bar{x}$, we also have in the state $\bar{z}_0$ that $\bar{\sigma}_{xy}^a = 0$ and there are no nonaffine flows, $\bar{v} = 0$ (Appendix B 2).
Combining Eqs. (C1)–(C3) with Eqs. (B1), (B2), and (B5) yields the following linearized dynamics:

$$ \bar{\partial}_t \delta s = -s'_0 \delta \bar{v}_x, \quad (C4) $$

$$ \begin{aligned} \bar{\partial}_t \delta c = & -c'_0 \delta \bar{v}_x + k_d(\delta s - \delta c) \\ & + D s_{ij}^{-2} \bar{\partial}_t \bar{\partial}_j \delta c - 2D(\bar{\partial}_x^2 c_0) \frac{\delta l_x}{l_{x0}^3}, \end{aligned} \quad (C5) $$

$$ \frac{d\delta l_x(t)}{dt} = \frac{2\operatorname{sgn}(\alpha)}{l_{x0}} \iint (\bar{\partial}_x^2 c_0) \delta c \, d\bar{x} d\bar{y} - \frac{dl_{x0}}{dt} \frac{\delta l_x}{l_{x0}}, \quad (C6) $$

where we have introduced the notation $s'_0 := \bar{\partial}_x s_0$ and $c'_0 := \bar{\partial}_x c_0$, and we used a partial integration in the linearized $l_x$ dynamics.

## 2. In codeforming Fourier space coordinates

In codeforming Fourier space, the linearized dynamics read

$$ \bar{\partial}_t \delta s = -s'_0 * \delta \bar{v}_x \quad (C7) $$

$$ \begin{aligned} \bar{\partial}_t \delta c = & -c'_0 * \delta \bar{v}_x + k_d(\delta s - \delta c) \\ & - D s_{ij}^{-2} \bar{q}_i \bar{q}_j \delta c + 2D \bar{q}_x^2 c_0 \frac{\delta l_x}{l_{x0}^3}, \end{aligned} \quad (C8) $$

$$ \frac{d\delta l_x(t)}{dt} = -\frac{2\operatorname{sgn}(\alpha)}{l_{x0}} \sum_{\bar{q}} \bar{q}_x^2 c_0^\dagger \delta c - \frac{dl_{x0}}{dt} \frac{\delta l_x}{l_{x0}}. \quad (C9) $$

Here $s'_0$ and $c'_0$ represent the codeforming Fourier transforms of $\bar{\partial}_x s_0$ and $\bar{\partial}_x c_0$, respectively. Furthermore, to rewrite the $\delta l_x$ dynamics, we have used Parseval's identity, where the dagger, $\dagger$, denotes the complex conjugate. We have moreover used the following definition of the convolution operator, *, in codeforming Fourier space:

$$ f(\bar{q}) * g(\bar{q}) := \sum_{\bar{q}'} f(\bar{q}') g(\bar{q} - \bar{q}'), \quad (C10) $$

where the sum is over all codeforming wave vectors $\bar{q}' = (\bar{q}'_x, \bar{q}'_y)$ with $\bar{q}'_i = 2\pi n'_i/L_i(0)$, where $n'_i \in \mathbb{Z}$ and $i \in \{x, y\}$. For a fixed system size, we ignore the $\delta l_x$ dynamics and the $\delta l_x$ term in the $\delta c$ dynamics.

To close Eqs. (C7)–(C9), we still need to insert the $x$ component of nonaffine flow field perturbation, $\delta \bar{v}_x$. For convenience, we will use *laboratory-frame* Fourier modes to compute $\delta \bar{v}_x$. Specifically, we use Eq. (B6), which expresses $\bar{v}$ in terms of the active stress tensor $\bar{\sigma}_{ij}^a$. We obtain to linear order

$$ \delta \bar{v}_x(q, \phi) = \frac{4i \sin \phi}{ql_{x0}} [\sin 2\phi \delta \bar{\sigma}_{xx}^a - \cos 2\phi \delta \bar{\sigma}_{xy}^a], \quad (C11) $$

where $q$ and $\phi$ denote the amplitude and orientation, respectively, of the laboratory-frame wave vector $q = q(\cos \phi, \sin \phi)$.

We further express $\delta \bar{\sigma}_{ij}^a$ in terms of $\delta c$ using Eq. (6). To linear order, in laboratory-frame Fourier space:

$$ \delta \bar{\sigma}_{xj}^a = \frac{i \operatorname{sgn}(\alpha)}{l_{x0}} c_0' * (q_j \delta c) \quad \text{for } j \in \{x, y\}. \quad (C12) $$

The other two components of the symmetric traceless tensor $\delta \bar{\sigma}_{ij}^a$ are obtained through the relations $\delta \bar{\sigma}_{yy}^a = -\delta \bar{\sigma}_{xx}^a$ and $\delta \bar{\sigma}_{yx}^a = \delta \bar{\sigma}_{xy}^a$. The convolution operator in laboratory-frame

Fourier coordinates in Eq. (C12) is defined analogously to Eq. (C10) as

$$ f(\boldsymbol{q}) * g(\boldsymbol{q}) := \sum_{\boldsymbol{q}'} f(\boldsymbol{q}') g(\boldsymbol{q} - \boldsymbol{q}'), \quad (C13) $$

where the sum is over all laboratory-frame wave vectors $\boldsymbol{q}' = (\boldsymbol{q}'_x, \boldsymbol{q}'_y)$ with $\boldsymbol{q}'_i = 2\pi n'_i/L_i(t)$, where $n'_i \in \mathbb{Z}$ and $i \in \{x, y\}$.

Combining (C11) and (C12), we obtain

$$ \delta \bar{v}_x = - \frac{4\operatorname{sgn}(\alpha) \sin \phi}{ql_{x0}^2} \{\sin 2\phi[c_0' * (q_x \delta c)] \\ - \cos 2\phi[c_0' * (q_y \delta c)]\}. \quad (C14) $$

The system of equations, Eqs. (C7)–(C9), and (C14), is now closed.

## 3. Numerical solution

In the special case of a linear source profile, we can solve the linearized dynamics analytically, because the gradients of $s_0$ and $c_0$ are constant, and thus the convolutions in Eqs. (C7), (C8), and (C14) collapse to a simple scaling factor. However, in general this is not the case, and we need to solve the linearized dynamics numerically.

To numerically compute the convolutions in Eqs. (C7), (C8), and (C14), we need to take into account that the numerical real-space representations of our fields are spatially discretized. To discuss the main ideas, we will discuss here the case of a single spatial dimension, $\bar{x}$, in codeforming coordinates. The generalization to laboratory frame coordinates and/or several dimensions will be straightforward. Specifically, we consider the periodic interval $\bar{x} \in [-L_x(0)/2, +L_x(0)/2)$ to be discretized by $N_x$ equal-length steps $\Delta \bar{x} = L_x(0)/N_x$, where $N_x$ is an even integer. Any function $f(\bar{x})$ is represented by a set of numbers $f_k$ with $k = -N_x/2, -N_x/2 + 1, \dots, 0, \dots, N_x/2 - 1$, such that $f_k = f(k \Delta \bar{x})$. As a consequence, the Fourier series are not infinite any more as in Eq. (A10), but instead cut off at some finite wave vectors. In particular, the numerical codeforming Fourier transform of $f$ is represented by a set of complex numbers $\tilde{f}_j$ with $j = -N_x/2, -N_x/2 + 1, \dots, 0, \dots, N_x/2 - 1$, where $\tilde{f}_j = f(\bar{q}_x = j\bar{q}_0)$ with $\bar{q}_0 := 2\pi/L_x(0)$. Eq. (A10) then becomes

$$ f_k = \sum_{j=-N_x/2}^{N_x/2-1} \tilde{f}_j e^{2\pi i (jk/N_x)}, \quad (C15) $$

where we used that $\bar{q}_0 \Delta \bar{x} = 2\pi/N_x$.

We define the convolution of two Fourier transforms, $\tilde{f}_j$ and $\tilde{g}_j$ as follows:

$$ (\tilde{f} * \tilde{g})_j = \sum_{j'=-N_x/2}^{N_x/2-1} \tilde{f}_{j'} \tilde{g}_{j-j'} \quad (C16) $$

for $j = -N_x/2, -N_x/2 + 1, \dots, 0, \dots, N_x/2 - 1$. Yet we see that the index $j - j'$ of $\tilde{g}$ can get outside of the range $-N_x/2, -N_x/2 + 1, \dots, N_x/2 - 1$. To prevent aliasing effects, we set $\tilde{g}_j \equiv 0$ whenever $j < -N_x/2$ or $j \ge N_x/2$. Note that this technique for dealiasing is akin to the zero-padding method [62,63], frequently used for spectral solvers.
To numerically solve the linearized dynamics, Eqs. (C7)–(C9), we use an exponential ansatz for the time dependence of the perturbations:

$$\delta s(\bar{q}, t) = \delta \hat{s}(\bar{q}) e^{\omega t}, \quad (C17)$$

$$\delta c(\bar{q}, t) = \delta \hat{c}(\bar{q}) e^{\omega t}, \quad (C18)$$

$$\delta l_x(t) = \delta \hat{l}_x e^{\omega t}. \quad (C19)$$

Inserting this in Eqs. (C7)–(C9) yields an eigenvalue problem:

$$\omega \delta \hat{z} = M \cdot \delta \hat{z}, \quad (C20)$$

where the components of $\delta \hat{z}$ are $\delta \hat{l}_x$ and the values of $\delta \hat{s}(\bar{q})$ and $\delta \hat{c}(\bar{q})$ for all $N_x N_y$ discrete values of the codeforming wave vector $\bar{q}$. Thus, $\delta \hat{z}$ has $2N_x N_y + 1$ components. The matrix $M$ usually couples not only $\delta \hat{s}$ and $\delta \hat{c}$ components for the same wave vector $\bar{q}$ to each other, but due to the convolutions in Eqs. (C7), (C8), and (C14), there is also mixing occurring across different wave vectors.

## APPENDIX D: FIXED-SIZE SYSTEM WITH LINEAR SOURCE

For linear source and signaling fields, the active flows, Eq. (C14), simplify to

$$\delta v_x = -\omega_{\text{act}} \delta c \quad (\text{D1})$$

with

$$\omega_{\text{act}} = 4 \operatorname{sgn}(\alpha) \sin^2 \phi, \quad (\text{D2})$$

where $\omega_{\text{act}}$ is the perturbation growth rate for the case of a boundary-provided gradient without source field [35]. With this, the linearized dynamics, Eqs. (C7) and (C8) become

$$\partial_t \delta s = \omega_{\text{act}} \delta c, \quad (\text{D3})$$

$$\partial_t \delta c = [\omega_{\text{act}} - Dq^2 - k_d] \delta c + k_d \delta s. \quad (\text{D4})$$

Thus, the convolutions in Eqs. (C7), (C8), and (C14) collapse to a constant prefactor. As a consequence, the solutions of the system are Fourier modes. In the following, we first discuss the limiting case of no diffusion before presenting the general solution.

### 1. No diffusion, $D = 0$

To solve the eigenvalue problem in this case, we define

$$\delta \Delta = \delta c - \delta s, \quad (\text{D5})$$

$$\delta c_{\text{eff}} = \delta c + \frac{k_d}{\omega_{\text{act}}} \delta s. \quad (\text{D6})$$

Inserting this into Eqs. (D3) and (D4), we obtain

$$\partial_t \delta \Delta = -k_d \delta \Delta, \quad (\text{D7})$$

$$\partial_t \delta c_{\text{eff}} = \omega_{\text{act}} \delta c_{\text{eff}}. \quad (\text{D8})$$

Hence, the difference between both signaling and source, $\delta \Delta$, decays at rate $k_d$, while a superposition of both, $\delta c_{\text{eff}}$,

grows with the same growth rate as the system without source [35].

### 2. General case

With an exponential Ansatz, where $\delta s \sim \delta c \sim e^{\omega t}$, we obtain an eigenvalue problem, whose solutions in $\omega$ are

$$\omega_{\pm} = \frac{1}{2} [T \pm \sqrt{T^2 + 4k_d \omega_{\text{act}}}], \quad (\text{D9})$$

where $T = \omega_{\text{act}} - Dq^2 - k_d$. The corresponding eigenvectors in $(\delta s, \delta c)$ are

$$v_+ = \begin{pmatrix} \omega_+ \\ \omega_{\text{act}} \end{pmatrix} \quad \text{and} \quad v_- = \begin{pmatrix} \omega_- \\ \omega_{\text{act}} \end{pmatrix}, \quad (\text{D10})$$

respectively.

## APPENDIX E: FIXED-SIZE SYSTEM WITH LOCALIZED SOURCE

### 1. Value of $s_b$ in dimensionless units

As we discussed in Sec. IIC, we set the value of $s_b$ such that the magnitude of the deformation rate of a system with $L_x = L_y = 1$ is one. This means that we set

$$s_b = \left[ \sum_{q_x} \frac{k_d^2 |m(q_x)|^2}{(k_d + Dq_x^2)^2} \right]^{-1/2}, \quad (\text{E1})$$

where $m(q_x)$ is the Fourier transform of the function

$$m(x) = \exp\left(\frac{\cos(q_0 x)}{(q_0 w)^2}\right), \quad (\text{E2})$$

which is defined on the interval $x \in [-1/2, 1/2]$, and where $q_0 = 2\pi/L_x$. Indeed, using Eqs. (7) and (22), one obtains that with the value of $s_b$ from Eq. (E1), a von Mises source profile, Eq. (21), and $L_x = L_y = 1$, the absolute value of the shear rate is one.

### 2. Convergence to a finite perturbation growth rate for $D \to \infty$

The large-diffusion limit can be understood similarly to Sec. III, where we assume that $Dq_{\text{min}}^2$ is faster than any other rate to eliminate the signaling perturbation dynamics. Specifically, only the source perturbation dynamics, Eq. (C7), remains, where the signaling perturbation $\delta c$ depends directly on the source perturbation through Eq. (19), which implies $\delta c \sim D^{-1}\delta s$. Moreover, for large $D$, Eq. (22) implies that $c_0 \sim D^{-1}s_0$, and thus $c'_0 \sim D^{-1}s'_0$.

In dimensionless units, we set $s_b$ such that the free deformation rate is one, which implies $s_0 \sim s_b \sim D$ [compare Eq. (E1)]. Thus, we have $c'_0 \sim D^0$. As a consequence, the perturbation of active stresses and active flows scale as $\delta \bar{v}_x \sim \delta \bar{\sigma}_{ij}^a \sim c'_0 \delta c \sim D^{-1}$ [compare Eqs. (C12) and (C14)]. The active flows enter in the advective term of the source field, Eq. (C7), where they are multiplied with the stationary source field, which scales as $s'_0 \sim s_b \sim D$. Hence, in the large-$D$ limit, the source dynamics does not depend on $D$ any more. As a consequence, the largest perturbation growth rate is constant (Fig. 6).

In dimensionful units, i.e., if we imposed a constant source magnitude $s_b \sim D^0$, we would have $c'_0 \sim D^{-1}$. Hence, the
<figure><img src="image_8.png" /><figcaption>FIG. 6. Largest perturbation growth rates $\omega_{\max}$ of a fixed-size system as a function of the signaling molecule diffusion. Parameter values, $k_d = 1$ and $L_x = L_y = 1$.</figcaption></figure>

perturbation growth rate would scale as $c'_0 \delta c \sim D^{-2}$, and the free deformation rate would scale as $\sim \bar{\sigma}_{xx}^a \sim (c'_0)^2 \sim D^{-2}$. Hence, the perturbation growth rate decreases with $D$ in the same way as the free deformation rate.

# APPENDIX F: DEFORMING SYSTEM WITH LOCALIZED SOURCE

## 1. System dimension dynamics

The dynamics of system dimension $L_x$, and thus box shear $l_x$, is given by Eq. (7).

We first discuss the case without diffusion, $D = 0$, where source and signaling fields, $s(\bar{x})$ and $c(\bar{x})$, are stationary in codeforming coordinates. In this case, the shear rate is essentially given by the integral over $(\partial_x c_0)^2$. Because the signaling profile, $c_0$, is stationary in codeforming coordinates, we have that for varying system size, $(\partial_x c_0)^2 \sim 1/l_x^2 (\bar{\partial}_x c_0)^2 \sim 1/l_x^2$. So, using dimensionless units (Sec. II C), we can write

$$ \frac{1}{l_x} \frac{dl_x}{dt} = - \frac{\operatorname{sgn}(\alpha)}{L_x^2(0)l_x^2}. \qquad (F1) $$

The solution of this is

$$ l_x(t) = \sqrt{1 - \frac{2\operatorname{sgn}(\alpha)t}{L_x(0)^2}}. \qquad (F2) $$

Thus, gradient-extensile systems ($\alpha < 0$) deform more and more slowly over time due to the widening of the signaling profile [green curves in Fig. 7(a)]. Meanwhile, gradient-contractile systems ($\alpha > 0$) deform faster and faster over time due to the compression of the signaling profile, until the deformation rate diverges at some critical time point $t_{crit} = L_x^2(0)/2$ [orange curves in Fig. 7(a)].

This box dimension dynamics is changed when allowing for finite diffusion, $D > 0$ [Fig. 7(b)]. Specifically, the shear rate of gradient-extensile systems first increases, before it

<figure><img src="image_9.png" /><figcaption>FIG. 7. System size dynamics of a freely deforming system in the limit $w \to \infty$, where the source has a cosine shape. Here $k_d = 1$. (a) Aspect ratios (left) and shear rates (right) over time for the gradient-extensile (green) and contractile (orange) systems, starting at initial aspect ratios of $L_x(0)/L_y(0) = 1/225$ and $L_x(0)/L_y(0) = 25$, respectively. Here the signaling molecule is nondiffusive, $D = 0$. (b) The same data as in panel (a) shown for a system with diffusive signaling molecule with $D = 10^{-2}$. Here the gradient-extensile (green) and gradient-contractile (orange) systems start from $L_x(0)/L_y(0) = 0.04$ and $L_x(0)/L_y(0) = 25$, respectively. Inset: Zoomed around the sharp drop in the gradient-contractile case.</figcaption></figure>

decreases like in the case without diffusion [compare green curves in Figs. 7(a) and 7(b), right panels]. In gradient contractile systems, the shear rate first increases, like in the case without diffusion before it rather abruptly drops [compare orange curves in Figs. 7(a) and 7(b), right panels and Fig. 7(b), inset].

To better understand these effects of diffusion, we first note that without any nonaffine flows (i.e., $\bar{v} = 0$), both signaling and source profiles keep their cosine profiles, $s(\mathbf{r}, t) = \hat{s}[\cos(2\pi x/L_x) + 1]$ and $c(\mathbf{r}, t) = \hat{c}(t) \cos(2\pi x/L_x) + c_{\text{offset}}(t)$, while only the signaling amplitude $\hat{c}(t)$ and offset $c_{\text{offset}}(t)$ vary over time. We can thus simplify the dynamics to

$$ \frac{d\hat{c}}{dt} = -k_d(\hat{c} - \hat{s}) - 4\pi^2 D \frac{\hat{c}}{L_x^2(0)l_x^2}, \qquad (F3) $$

$$ \frac{1}{l_x} \frac{dl_x}{dt} = - \frac{2\pi^2 \operatorname{sgn}(\alpha) \hat{c}^2}{L_x^2(0)l_x^2}, \qquad (F4) $$

where $\hat{s}$ is constant since $\bar{v} = 0$, and the offset $c_{\text{offset}}(t)$ does not matter for the box dimension dynamics. Equation (F4) implies that the shear rate scales as $\sim \hat{c}^2/l_x^2$. Further, the diffusion term in Eq. (F3) can be neglected as long as $l_x \gtrsim \sqrt{D/k_d/L_x(0)}$, and consequentially $\hat{c} \simeq \hat{s}$. In other words, the box shear rate for large $l_x$ is the same as in the case without diffusion, where the shear rate scales as $\sim l_x^{-2}$, for both gradient-extensile and gradient-contractile systems, consistent with Figs. 7(a) and 7(b). Conversely, taking Eq. (F3) to the adiabatic limit of large diffusion $D$ and/or small box shear $l_x$, we have $\hat{c} \sim (k_d/D)L_x^2(0)l_x^2 \hat{s} \sim l_x^2$. Thus, in the limit of small box shear $l_x$, the shear rate scales as $\sim \hat{c}^2/l_x^2 \sim l_x^2$, i.e., it increases with $l_x$, again consistent with Figs. 7(a) and 7(b).
## 2. Affine diffusion-shear instability

To discuss the affine diffusion-shear instability, we start from the linearized dynamics in Fourier space, Eqs. (C7)–(C9), where we set nonaffine flows to zero, $\delta\bar{v} = 0$, we ignore any variations of the source field, $\delta s = 0$, and we ignore any dependence of the perturbations. We obtain from Eqs. (C8) and (C9) with an exponential Ansatz for the time dependency, $\delta c(\bar{q}_x, t) = \delta\bar{c}(\bar{q}_x)e^{\omega t}$ and $\delta l_x(t) = \delta\bar{l}_x e^{\omega t}$:

$$ 
\begin{aligned}
\omega\delta\bar{c}(\bar{q}_x) &= -k_d(1 + \lambda^2 q_x^2)\delta\bar{c}(\bar{q}_x) + 2k_d\lambda^2 q_x^2 c_0(\bar{q}_x)\frac{\delta\bar{l}_x}{l_x}, \quad (\text{F5}) \\
\omega\frac{\delta\bar{l}_x}{l_x} &= -2\operatorname{sgn}(\alpha)\sum_{\bar{q}_x} q_x^2 c_0^\dagger(\bar{q}_x)\delta\bar{c}(\bar{q}_x) \\
&\quad + \operatorname{sgn}(\alpha)\frac{\delta\bar{l}_x}{l_x}\sum_{\bar{q}_x} q_x^2 |c_0(\bar{q}_x)|^2.
\end{aligned}
\quad (\text{F6})
$$

Here we also substituted the diffusion length scale $\lambda = \sqrt{D/k_d}$ and the laboratory-frame wave vector $q_x = \bar{q}_x/l_x$. If and only if this linear system of equations has a positive eigenvalue $\omega$, then there is an instability.

### a. Cosine source profile

For the case of a cosine profile, the only nonzero values of $c_0(\bar{q}_x)$ are $c_0 \equiv c_0(\pm\bar{q}_0)$ with $\bar{q}_0 = 2\pi/L_x(0)$. Then Eqs. (F5) and (F6) become

$$ \omega\delta\bar{c} = -k_d(1 + \lambda^2 q_0^2)\delta\bar{c} + 2k_d\lambda^2 q_0^2 c_0\frac{\delta\bar{l}_x}{l_x}, \quad (\text{F7}) $$

$$ \omega \frac{\delta \bar{l}_x}{l_x} = -4\operatorname{sgn}(\alpha) q_0^2 c_0^\dagger \delta \bar{c} + 2\operatorname{sgn}(\alpha) q_0^2 |c_0|^2 \frac{\delta \bar{l}_x}{l_x}, \quad (\text{F8}) $$

where $q_0(t) := \bar{q}_0/l_x(t) = 2\pi/L_x(t)$. There are two eigenvalues, $\omega_1$ and $\omega_2$ with $\omega_1 \le \omega_2$, and the largest is given by

$$ \omega_2 = \frac{1}{2}[T + \sqrt{T^2 - 4D}], \quad (\text{F9}) $$

where $T := -k_d(1 + \lambda^2 q_0^2) + 2\operatorname{sgn}(\alpha)q_0^2|c_0|^2$ and $D := -2\operatorname{sgn}(\alpha)k_d q_0^2|c_0|^2(1 - 3\lambda^2 q_0^2)$ are trace and determinant, respectively, of the matrix on the right-hand side of Eqs. (F7) and (F8).

In the gradient-extensile case, $\operatorname{sgn}(\alpha) = -1$, we obtain that $\omega_1 + \omega_2 = T < 0$. Thus, at least one eigenvalue is negative. In addition, the sign of the product, $\omega_1\omega_2 = D$, will indicate whether both eigenvalues are negative or at least one is positive. The product is given by $\omega_1\omega_2 = D = 2k_d q_0^2|c_0|^2(1 - 3\lambda^2 q_0^2)$. Hence, there is a positive eigenvalue if and only if $\omega_1\omega_2 < 0$, which is equivalent to $3\lambda^2 q_0^2 > 1$, and thus to $2\pi\lambda\sqrt{3} > L_x$.

### b. Arbitrary source profile

For the case of an arbitrary source profile $s_0$, we focus on the gradient-extensile case, $\operatorname{sgn}(\alpha) = -1$. We discuss the existence of positive solutions $\omega$ by first expressing $\delta\bar{c}(\bar{q}_x)$ in terms of $\delta\bar{l}_x$ and $\omega$ using Eq. (F5):

$$ \delta\bar{c}(\bar{q}_x) = \frac{2\lambda^2 q_x^2 c_0(\bar{q}_x)}{1 + \omega/k_d + \lambda^2 q_x^2} \frac{\delta\bar{l}_x}{l_x}. \quad (\text{F10}) $$

Insertion into Eq. (F6) yields

$$ 
\begin{aligned}
\omega &= \sum_{\bar{q}_x} q_x^2 |c_0(\bar{q}_x)|^2 \frac{4\lambda^2 q_x^2}{1 + \omega/k_d + \lambda^2 q_x^2} - \sum_{\bar{q}_x} q_x^2 |c_0(\bar{q}_x)|^2 \\
&= \sum_{\bar{q}_x} q_x^2 |c_0(\bar{q}_x)|^2 \frac{3\lambda^2 q_x^2 - 1 - \omega/k_d}{1 + \omega/k_d + \lambda^2 q_x^2}.
\end{aligned}
\quad (\text{F11})
$$

To solve this self-consistent equation for $\omega$, we look for the zeros of the function

$$ f_{\lambda}(\omega) := -\omega + \sum_{\bar{q}_x} q_x^2 |c_0(\bar{q}_x)|^2 \frac{3\lambda^2 q_x^2 - 1 - \omega/k_d}{1 + \omega/k_d + \lambda^2 q_x^2}. \quad (\text{F12}) $$

Specifically, we now show that iff $f_{\lambda}(0) > 0$, there exists a positive solution in $\omega$, i.e., an $\omega > 0$ for which $f_{\lambda}(\omega) = 0$. This means that there is an affine diffusion-shear instability iff $f_{\lambda}(0) > 0$. To show this, we first assume that there is an $\omega > 0$ for which $f_{\lambda}(\omega) = 0$. Then we have

$$ 
\begin{aligned}
0 < \omega &= \sum_{\bar{q}_x} q_x^2 |c_0(\bar{q}_x)|^2 \frac{3\lambda^2 q_x^2 - 1 - \omega/k_d}{1 + \omega/k_d + \lambda^2 q_x^2} \\
&< \sum_{\bar{q}_x} q_x^2 |c_0(\bar{q}_x)|^2 \frac{3\lambda^2 q_x^2 - 1}{1 + \lambda^2 q_x^2} \\
&= f_{\lambda}(0).
\end{aligned}
\quad (\text{F13})
$$

In these transformations, we have subsequently used $\omega > 0$, $f_{\lambda}(\omega) = 0$, $\omega > 0$, and the definition of $f_{\lambda}(\omega)$. Thus, if there is an $\omega > 0$ for which $f_{\lambda}(\omega) = 0$, then $f_{\lambda}(0) > 0$.

Second, we assume that $f_{\lambda}(0) > 0$. To show that there is a positive zero of $f_{\lambda}$, we need to show that there is a positive $\hat{\omega} > 0$ for which $f_{\lambda}(\hat{\omega}) < 0$. Then, the Intermediate Value Theorem implies that there is an $\omega$ with $0 < \omega < \hat{\omega}$ such that $f_{\lambda}(\omega) = 0$. Indeed, for $\hat{\omega} = 3 \sum_{\bar{q}_x} q_x^2 |c_0(\bar{q}_x)|^2$, we find that

$$ 
\begin{aligned}
f_{\lambda}(\hat{\omega}) &= -3 \sum_{\bar{q}_x} q_x^2 |c_0(\bar{q}_x)|^2 \\
&\quad + \sum_{\bar{q}_x} q_x^2 |c_0(\bar{q}_x)|^2 \frac{3\lambda^2 q_x^2 - 1 - \hat{\omega}/k_d}{1 + \hat{\omega}/k_d + \lambda^2 q_x^2} \\
&= \sum_{\bar{q}_x} q_x^2 |c_0(\bar{q}_x)|^2 \frac{-2(1 + \hat{\omega}/k_d)}{1 + \hat{\omega}/k_d + \lambda^2 q_x^2} \\
&< 0.
\end{aligned}
\quad (\text{F14})
$$

Thus, if $f_{\lambda}(0) > 0$ then $f_{\lambda}$ has a positive zero and the affine diffusion-shear instability exists. Taken together, the sign of $f_{\lambda}(0)$ directly indicates the existence of an instability.

Finally, we discuss the sign of $f_{\lambda}(0)$. Using Eq. (22), we can express $f_{\lambda}(0)$ in terms of the source profile:

$$ f_{\lambda}(0) = \frac{1}{\lambda^2} \sum_{\bar{q}_x} |s_0(\bar{q}_x)|^2 \frac{\lambda^2 q_x^2 (3\lambda^2 q_x^2 - 1)}{(1 + \lambda^2 q_x^2)^3}. \quad (\text{F15}) $$

Choosing for $s_0$ the von Mises profile in Eq. (21) with finite width $w$, we numerically obtain $f_{\lambda}(0) > 0$ iff $\lambda > \text{gw}$ with $g \approx 0.639$ in the limit where $w \ll L_x$ (i.e., where the von Mises profile becomes a Gaussian with standard deviation $w$).
However, the precise value of g depends on the shape of the profile $s_0$. For a rectangular profile of $s_0$, we numerically even find that the instability always occurs.

## APPENDIX G: NONLINEARITY IN THE ACTIVE STRESS

We consider the system of equations, Eqs. (1)–(6), where we introduce a nonlinearity into the active stress, Eq. (6):

$$ \bar{\sigma}_{ij}^a = \alpha h([\partial_k c]^2) [(\partial_i c)(\partial_j c) - \frac{1}{2} (\partial_k c)^2 \delta_{ij}], \quad (G1) $$

where $z \mapsto h(z)$ is some real function expressing the nonlinearity. In the following, we assume that $h(z) > 0$ for any real $z$, and that the function $z \mapsto zh(z)$ is monotonously increasing, i.e., the magnitude of the active stress increases monotonously with the magnitude of scalar field gradient.

Then the dimensionless, linearized dynamics become

$$ \bar{\partial}_t \delta s = -s'_0 \delta \bar{v}_x, \quad (G2) $$

$$ \begin{aligned} \bar{\partial}_t \delta c = & -c'_0 \delta \bar{v}_x + k_d (\delta s - \delta c) \\ & + D s_{ij}^{-2} \bar{\partial}_t \bar{\partial}_j \delta c - 2D (\bar{\partial}_x^2 c_0) \frac{\delta l_x}{l_{x0}^3}, \end{aligned} \quad (G3) $$

$$ \begin{aligned} \frac{d\delta l_x(t)}{dt} = & -2\operatorname{sgn}(\alpha) \iint [h(z) + zh'(z)] \frac{c_0'}{l_{x0}} \delta c' \, d\bar{x}d\bar{y} \\ & + \operatorname{sgn}(\alpha)(I_0 + 2\Delta I) \delta l_x, \end{aligned} \quad (G4) $$

where $\delta c' := \bar{\partial}_x \delta c$, and

$$ I_0 := \iint zh(z) d\bar{x}d\bar{y}, \quad (G5) $$

$$ \Delta I := \iint z^2 h'(z) d\bar{x}d\bar{y} \quad (G6) $$

with $z(\bar{r}) := (c'_0(\bar{x})/l_{x0})^2$.

We further focus on the gradient-extensile case, sgn($\alpha$) = -1, without nonaffine flows, $\bar{v} = 0$, and, for simplicity, we consider the limit $w \to \infty$, such that $s_0(\bar{x}) \sim \cos(\bar{q}_0 \bar{x})$, and thus also $c_0(\bar{x}) = \hat{c}_0 \cos(\bar{q}_0 \bar{x})$ with $\bar{q}_0 = 2\pi/L_x(0)$ and some constant prefactor $\hat{c}_0$. As a consequence, based on Eq. (G3), all Fourier modes of $\delta c$ decay except for $\delta c = \delta \hat{c} \cos(\bar{q}_0 \bar{x})$. Hence, $\delta c(\bar{x}) = (\delta \hat{c}/\hat{c}_0)c_0(\bar{x})$. As a consequence, Eqs. (G3) and (G4) simplify to

$$ \frac{d\delta \hat{c}}{dt} = -(k_d + Dq_0^2)\delta \hat{c} + 2Dq_0^2 \hat{c}_0 \frac{\delta l_x}{l_{x0}}, \quad (G7) $$

$$ \frac{d\delta l_x(t)}{dt} = 2(I_0 + \Delta I)l_{x0} \frac{\delta \hat{c}}{\hat{c}_0} - (I_0 + 2\Delta I)\delta l_x. \quad (G8) $$

This is a linear system of differential equations with constant coefficients. It is linearly stable if both eigenvalues, $\omega_1$ and $\omega_2$, of the coefficient matrix are negative. This implies that the trace of the matrix needs to be negative: $0 > \omega_1 + \omega_2 = -k_d - Dq_0^2 - I_0 - 2\Delta I$, and that the determinant of the matrix needs to be positive: $0 < \omega_1\omega_2 = (k_d + Dq_0^2)(I_0 + 2\Delta I) - 4Dq_0^2(I_0 + \Delta I) = k_d(I_0 + 2\Delta I) - Dq_0^2(3I_0 + 2\Delta I)$. For the system to be stable also for large diffusion, i.e., for $k_d \ll Dq_0^2$, the latter condition implies $-3I_0 - 2\Delta I > 0$. However, we have assumed that $h(z) > 0$, implying $I_0 > 0$, and that $z \mapsto zh(z)$ increases monotonically, implying that $h(z) + zh'(z) \ge 0$ and thus $I_0 + \Delta I \ge 0$. This contradicts determinant condition for a stable system. Taken together, this means that the diffusion-shear instability arises for large diffusion for any nonlinearity $h$ in the active stress. We expect this result to generalize also to finite source widths $w$; higher-order Fourier modes in $c_0$ would increase prefactor of the first term in Eq. (G8) and thus make the system more unstable.

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

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

[3] B. Bénazéraf, P. François, R. E. Baker, N. Denans, C. D. Little, and O. Pourquié, A random cell motility gradient downstream of FGF controls elongation of an amniote embryo, *Nature* (London) **466**, 248 (2010).

[4] A. Shindo, Models of convergent extension during morphogenesis, *WIREs Dev. Biol.* **7**, e293 (2018).

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

[6] C. M. Karner, R. Chirumamilla, S. Aoki, P. Igarashi, J. B. Wallingford, and T. J. Carroll, Wnt9b signaling regulates planar cell polarity and kidney tubule morphogenesis, *Nat. Genet.* **41**, 793 (2009).

[7] A. Saxena, B. Denholm, S. Bunt, M. Bischoff, K. VijayRaghavan, and H. Skaer, Epidermal growth factor signaling controls myosin II planar polarity to orchestrate convergent extension movements during *Drosophila* tubulogenesis, *PLoS Biol.* **12**, e1002013 (2014).

[8] R. Etournay, M. Popović, M. Merkel, A. Nandi, C. Blasse, B. Aigouy, H. Brandl, G. Myers, G. Salbreux, F. Jülicher, and S. Eaton, Interplay of cell dynamics and epithelial tension during morphogenesis of the *Drosophila* pupal wing, *eLife* **4**, e07090 (2015).

[9] S. Hopyan, Biophysical regulation of early limb bud morphogenesis, *Dev. Biol.* **429**, 429 (2017).

[10] H. Tao, M. Zhu, K. Lau, O. K. W. Whitley, M. Samani, X. Xiao, X. X. Chen, N. A. Hahn, W. Liu, M. Valencia, M. Wu, X. Wang, K. D. Fenelon, C. C. Pasiliao, D. Hu, J. Wu, S. Spring, J. Ferguson, E. P. Karuna, R. M. Henkelman et al., Oscillatory cortical forces promote three dimensional cell intercalations that shape the murine mandibular arch, *Nat. Commun.* **10**, 1703 (2019).

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

[12] F. Bosveld, I. Bonnet, B. Guirao, S. Tlili, Z. Wang, A. Petitalot, R. Marchand, P.-L. Bardet, P. Marcq, F. Graner, and Y. Bellaiche, Mechanical control of morphogenesis by
fat/dachsous/four-jointed planar cell polarity pathway, *Science* **336**, 724 (2012).

[13] C. Collinet, M. Rauzi, P. F. Lenne, and T. Lecuit, Local and tissue-scale forces drive oriented junction growth during tissue extension, *Nat. Cell Biol.* **17**, 1247 (2015).

[14] M. Behrndt, G. Salbreux, P. Campinho, R. Hauschild, F. Oswald, J. Roensch, S. W. Grill, and C. P. Heisenberg, Forces driving epithelial spreading in zebrafish gastrulation, *Science* **338**, 257 (2012).

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

[16] A. Stokkermans, A. Chakrabarti, K. Subramanian, L. Wang, S. Yin, P. Moghe, P. Steenbergen, G. Mönke, T. Hiiragi, R. Prevedel, L. Mahadevan, and A. Ikmi, Muscular hydraulics drive larva-polyp morphogenesis, *Curr. Biol.* **32**, 4707 (2022).

[17] E. W. Gehrels, B. Chakraborty, M.-E. Perrin, M. Merkel, and T. Lecuit, Curvature gradient drives polarized tissue flow in the *Drosophila* embryo, *Proc. Natl. Acad. Sci. USA* **120**, e2214205120 (2023).

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

[19] M. Serra, G. S. Nájera, M. Chuai, A. M. Plum, S. Santhosh, V. Spandan, C. J. Weijer, and L. Mahadevan, A mechanochemical model recapitulates distinct vertebrate gastrulation modes, *Sci. Adv.* **9**, eadh8152 (2023).

[20] S. Gsell, S. Tlili, M. Merkel, and P.-F. Lenne, Marangoni-like tissue flows enhance symmetry breaking of embryonic organoids, *Nat. Phys.* **21**, 644 (2025).

[21] K. Barrett, S. Anand, V. Thome, L. Kodjabachian, M. Merkel, and P.-F. Lenne, Epithelial-mesenchymal interface guides cell shapes and axis elongation in embryonic explants, *Development* **152**, dev204344 (2025).

[22] N. H. Claussen, F. Brauns, and B. I. Shraiman, A geometric-tension-dynamics model of epithelial convergent extension, *Proc. Natl. Acad. Sci. USA* **121**, e2321928121 (2024).

[23] A. Ioratim-Uba, T. B. Liverpool, and S. Henkes, Mechanochemical active feedback generates convergence extension in epithelial tissue, *Phys. Rev. Lett.* **131**, 238301 (2023).

[24] R. A. Simha and S. Ramaswamy, Hydrodynamic fluctuations and instabilities in ordered suspensions of self-propelled particles, *Phys. Rev. Lett.* **89**, 058101 (2002).

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

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

[27] S. von der Hardt, J. Bakkers, A. Inbal, L. Carvalho, L. Solnic-Krezel, C. P. Heisenberg, and M. Hammerschmidt, The bmp gradient of the zebrafish gastrula guides migrating lateral cells by regulating cell-cell adhesion, *Curr. Biol.* **17**, 475 (2007).

[28] K. W. Rogers and A. F. Schier, Morphogen gradients: From generation to interpretation, *Annu. Rev. Cell Dev. Biol.* **27**, 377 (2011).

[29] P. Müller, K. W. Rogers, S. R. Yu, M. Brand, and A. F. Schier, Morphogen transport, *Development* **140**, 1621 (2013).

[30] G. Tkačik, J. O. Dubuis, M. D. Petkova, and T. Gregor, Positional information, positional error, and readout precision in morphogenesis: A mathematical framework, *Genetics* **199**, 39 (2015).

[31] L. S. Mosby, A. E. Bowen, and Z. Hadjivasiliou, Morphogens in the evolution of size, shape and patterning, *Development* **151**, dev202412 (2024).

[32] H. Ninomiya, R. P. Elinson, and R. Winklbauer, Anteroposterior tissue polarity links mesoderm convergent extension to axial patterning, *Nature (London)* **430**, 364 (2004).

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

[34] Z. Wang, M. C. Marchetti, and F. Brauns, Patterning of morphogenetic anisotropy fields, *Proc. Natl. Acad. Sci. USA* **120**, e2220167120 (2023).

[35] M. Ibrahimi and M. Merkel, Deforming polar active matter in a scalar field gradient, *New J. Phys.* **25**, 013022 (2023).

[36] M. F. Lefebvre, N. H. Claussen, N. P. Mitchell, H. J. Gustafson, and S. J. Streichan, Geometric control of myosin-II orientation during axis elongation, *eLife* **12**, e78787 (2023).

[37] A. M. Plum and M. Serra, Morphogen patterning in dynamic tissues, *PRX Life* **3**, 043009 (2025).

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

[39] R. Mateus, L. Holtzer, C. Seum, Z. Hadjivasiliou, M. Dubois, F. Jülicher, and M. Gonzalez-Gaitán, BMP signaling gradient scaling in the zebrafish pectoral fin, *Cell Rep.* **30**, 4292 (2020).

[40] M. Romanova-Michaelides, Z. Hadjivasiliou, D. Aguilar-Hidalgo, D. Basagiannis, C. Seum, M. Dubois, F. Jülicher, and M. Gonzalez-Gaitán, Morphogen gradient scaling by recycling of intracellular dpp, *Nature (London)* **602**, 287 (2022).

[41] O. Wartlick, A. Kicheva, and M. González-Gaitán, Morphogen gradient formation, *Cold Spring Harbor Perspect. Biol.* **1**, a001255 (2009).

[42] Y. Wang, X. Wang, T. Wohland, and K. Sampath, Extracellular interactions and ligand degradation shape the nodal morphogen gradient, *eLife* **5**, e13879 (2016).

[43] A. Huang, C. Amourda, S. Zhang, N. S. Tolwinski, and T. E. Saunders, Decoding temporal interpretation of the morphogen bicoid in the early *Drosophila* embryo, *eLife* **6**, 1 (2017).

[44] A. C. Paré, A. Vichas, C. T. Fincher, Z. Mirman, D. L. Farrell, A. Mainieri, and J. A. Zallen, A positional toll receptor code directs convergent extension in *Drosophila*, *Nature (London)* **515**, 523 (2014).

[45] M. A. Benton, M. Pechmann, N. Frey, D. Stappert, K. H. Conrads, Y. T. Chen, E. Stamataki, A. Pavlopoulos, and S. Roth, Toll genes have an ancestral role in axis elongation, *Curr. Biol.* **26**, 1609 (2016).

[46] D. Jia, Q. Xu, Q. Xie, W. Mio, and W. M. Deng, Automatic stage identification of *Drosophila* egg chamber based on DAPI images, *Sci. Rep.* **6**, 18850 (2016).

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

[48] D. D. Iwaki, K. A. Johansen, J. B. Singer, and J. A. Lengyel, Drumstick, bowl, and lines are required for patterning and cell
rearrangement in the *Drosophila* embryonic hindgut, *Dev. Biol.* **240**, 611 (2001).

[49] M. L. Williams and L. Solnica-Krezel, Nodal and planar cell polarity signaling cooperate to regulate zebrafish convergence and extension gastrulation movements, *eLife* **9**, e54445 (2020).

[50] B. Steventon, F. Duarte, R. Lagadec, S. Mazan, J. F. Nicolas, and E. Hirsinger, Species-specific contribution of volumetric growth and tissue convergence to posterior body elongation in vertebrates, *Development (Cambridge)* **143**, 1732 (2016).

[51] A. Kicheva, P. Pantazis, T. Bollenbach, Y. Kalaidzidis, T. Bittig, F. Jülicher, and M. González-Gaitán, Kinetics of morphogen gradient formation, *Science* **315**, 521 (2007).

[52] A. Tiribocchi, R. Wittkowski, D. Marenduzzo, and M. E. Cates, Active model H: Scalar active matter in a momentum-conserving fluid, *Phys. Rev. Lett.* **115**, 188302 (2015).

[53] T. R. Kirkpatrick and J. K. Bhattachacherjee, Driven active matter: Fluctuations and a hydrodynamic instability, *Phys. Rev. Fluids* **4**, 024306 (2019).

[54] B. Gao, Wnt regulation of planar cell polarity (PCP), *Curr. Topics Dev. Biol.* **101**, 263 (2012).

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

[56] See Supplemental Material at http://link.aps.org/supplemental/10.1103/yq1s-p2z1 for movies showing the simulations.

[57] S. S. Lienkamp, K. Liu, C. M. Karner, T. J. Carroll, O. Ronneberger, J. B. Wallingford, and G. Walz, Vertebrate kidney tubules elongate using a planar cell polarity-dependent, rosette-based mechanism of convergent extension, *Nat. Genet.* **44**, 1382 (2012).

[58] O. Ossipova, K. Kim, and S. Y. Soko, Planar polarization of Vangl2 in the vertebrate neural plate is controlled by Wnt and Myosin II signaling, *Biology Open* **4**, 722 (2015).

[59] M. Ibrahimi, Robustness of active anisotropic deformation in developing biological tissues, Ph.D. thesis, Aix-Marseille Université, 2022.

[60] M. Lefebvre, J. Colen, N. Claussen, F. Brauns, M. Raich, N. Mitchell, M. Fruchart, V. Vitelli, and S. J. Streichan, Learning a conserved mechanism for early neuroectoderm morphogenesis, *arXiv:2405.18382v1*.

[61] M. Ibrahimi and M. Merkel, Scripts for “Stabilization of active tissue deformation by a dynamic signaling gradient”, *PRX Life* (2025), doi:10.24433/CO.8370709.v1.

[62] M. Roberts and J. C. Bowman, Dealiased convolutions for pseudospectral simulations, *J. Phys.: Conf. Ser.* **318**, 072037 (2011).

[63] J. C. Bowman and M. Roberts, Efficient dealiased convolutions without padding, *SIAM J. Sci. Comput.* **33**, 386 (2011).