# Passive defect driven morphogenesis in nematic membranes

D. J. G. Pearce, C. Thibault, Q. Chaboche, C. Blanch-Mercader

► To cite this version:

D. J. G. Pearce, C. Thibault, Q. Chaboche, C. Blanch-Mercader. Passive defect driven morphogenesis in nematic membranes. Physical Review Letters, 2024, 134 (1), pp.018402. 10.1103/Phys-RevLett.134.018402. hal-04780983v2

HAL Id: hal-04780983

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

Submitted on 26 Mar 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
Passive defect driven morphogenesis in nematic membranes

D. J. G. Pearce,¹,* C. Thibault,² Q. Chaboche,² and C. Blanch-Mercader²,†

¹Dept. of Theoretical Physics, University of Geneva, 1211 Geneva, Switzerland

²Laboratoire Physico-Chimie Curie, Institut Curie, Université PSL, Sorbonne Université, CNRS UMR168, F-75248 Paris, France

Orientational order is a common feature of many biological and synthetic materials. Topological defects are discontinuities in this order that are often coupled to geometric features of the materials. Here, we study the equilibrium shapes of fluid membranes featuring a +1 topological defect as a model for morphogenesis. We show, through simulation and analytic calculation, that the membrane can spontaneously deform toward a conical shape with a defect at its apex. We show that the relative stability of the deformation is controlled by the balance of the elastic parameters. When boundary constraints are introduced, we observe three distinct modes of deformation. These deformation modes take advantage of the way in which splay, twist and bend distortions of the director field can be exchanged on a curved surface. Finally we demonstrate inverted solutions. Our findings demonstrate a mechanism for passive defect driven morphogenesis as well as the fusion of +1/2 topological defect pairs on deformable surfaces.

The study of shape shifting materials is often inspired by biological systems and has a broad range of applications. For example, deformable solids that have a prescribed orientational field, such as elastomers [1–3] or inflatable structures [4, 5], can achieve families of morphologies by differential growth [6–10] and have applications in soft robotics[11]. Here, we focus on fluid membranes with nematic order. Unlike the previous cases, these membranes can, in addition, flow, remodel and adapt both their orientational field and shape.

Fluid membranes achieve their shapes by minimizing a free-energy subjected to physical constraints [12–15]. When orientational order is included, the stresses generated by topological defects can render a flat surface unstable to out-of-plane perturbations [16]. In addition, the geometrical properties of the surface, such as the extrinsic curvature, can in turn influence the dynamics of the orientational field and its equilibrium configurations [17–24]. Indeed, for prescribed geometries, couplings between the nematic field and the intrinsic or the extrinsic geometry have been shown to induce symmetry-breaking of orientational configurations, or influence topological defect dynamics, [17–19, 25–27]. Recent studies extended these works by, for instance, including out-of-equilibrium processes, such as active stresses or active torques, or varying surface geometry and topology [28–36].

We study rotationally symmetric systems, which therefore feature a +1 topological defect at their center; +1 topological defects have been associated with geometrical changes in natural and synthetic systems [37–40]. We use cylindrical coordinates to define these shapes, where r is the radial coordinate, θ is the azimuthal coordinate and ζ is the axial coordinate, Fig. 1a. The outer circular boundary is placed at r = R; without loss of generality, we set the vertical offset by ζ(R) = 0 and R = 1. The orientation of the nematic field is described by a director field ˆn that represents the averaged local orientation on the surface. We consider that the director field is tangential

to the surface. In addition, we consider that the system is deep into the nematic phase and impose | ˆn | = 1. This allows the director field to be defined by a scalar phase ψ(r), which corresponds to the angle between the director field and the curvilinear radial direction. The cases ψ = 0 corresponds to the aster, 0 < ψ < π/2 to a spiral, and ψ = π/2 to the vortex.

The two-dimensional free-energy of a fluid membrane with nematic order is given by

$$
\mathcal{F} = \int_{\mathcal{M}} \left\{ k_B \mathcal{H}^2 + \sigma + k_1 (\nabla \cdot \hat{n})^2 + k_2 (\hat{n} \cdot (\nabla \times \hat{n}))^2 + k_3 (\hat{n} \times (\nabla \times \hat{n}))^2 \right\} da. \quad (1)
$$

The first term is the bending energy with mean curvature H, and the second term represents surface tension. We disregard anisotropies in bending energy dependent on ˆn. The other terms are the Frank free-energy associated respectively with splay, twist and bend distortions of the director field ˆn [41, 42]. The corresponding elastic coefficients are: the bending rigidity k_B, the surface tension σ, and the reduced Frank constants k_1, k_2, and k_3 that are proportional to the membrane thickness. Because the director field is parallel to the surface, the effects of the saddle-splay distortions with elastic constant k_24 can be absorbed in a redefinition of the other Frank constants, [18]. As shown in Refs. [17, 18] and subsequently extended for more general cases in [43], the thin-film limit of the Frank free-energy results in contributions that couple the director field with both the intrinsic and the extrinsic geometry. In our case, Eq. (1) takes the form derived in Refs. [17, 18] and the expression of the free-energy (1) for the special case of a surface of revolution is derived in Sec. 1.1 [44]. Other descriptions for fluid surfaces with orientational order neglected the couplings with the extrinsic geometry [16, 25, 26]. For a discussion on the thin-film approximation of liquid crystals, we refer to Ref. [43].
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Geometrical setup. The red line shows the height profile, $\zeta(r)$, and black dashes show the orientation of the director, $\psi(r)$. (b) The height of a membrane featuring a single aster defect with equal Frank constants $k_1=k_2=k_3=k/3$ as a function of $k$, surface tension $\sigma$ and bending rigidity $k_B$ with the constraint $\sigma + k_B + k = 1$. The magenta line corresponds to Eq. 4. From this point forward we fix $\sigma = k_B = 1/10$ and $k_1 + k_2 + k_3 = 1$ in all results unless otherwise stated. (c) Height of a membrane as a function of the prescribed uniform phase, $\psi$, and ratio between the splay and bend coefficients, ($k_2 = 0$). The green line corresponds to the predicted transition to conical surfaces, see SI. (d)&(e) Examples of equilibrium surfaces with parameters given by: white star ($k_1 = k_3$, $\psi = 0$), red star ($k_3 = 1$, $\psi = \pi/8$). (d) An aster on a conical surface, and (e) a spiral on a concave surface. Blue (red) regions: negative (positive) Gaussian curvature.</figcaption></figure>

We combine analytical and numerical approaches to study equilibrium configurations of the free-energy (1). The former approach is restricted to director fields with a uniform phase $\psi$, see Sec. 1 and 2 in [44]. The latter approach discretizes the functions $\zeta(r)$ and $\psi(r)$, Fig. 1a, and uses a Monte-Carlo algorithm, see Sec. 3 [44].

We first consider the case where all the reduced Frank constants are equal, i.e. the one-constant approximation $k_1=k_2=k_3=k/3$. The energy scale is set by the constraint $k+k_B+\sigma=1$. Fig. 1b shows the height at the center $\zeta(0)$ of the equilibrium states that were found for varying elastic parameters. When $k_B$ or $\sigma$ dominate, the minimal states are flat discs which minimize surface area and curvature. Because the Frank constants are equal, $k_1=k_3$, all director configurations with a constant phase $\psi$ have equal energy, and thus are minimal states, [42]. However, when $k$ dominates, the minimal state changes to an approximately conical surface with phase $\psi = 0$, Fig. 1d. The morphological transition occurs via a symmetry-breaking process during which the director adopts an aster configuration.

To better understand this spontaneous transition from flat to conical surfaces, we derived analytically the normal force balance equation for a surface of revolution

with an embedded director field with uniform phase $\psi$, see Sec. 1.2 in [44] (see [27] for a general case). In the special case $\sigma=0$, this nonlinear ODE has a set of exact non-trivial solutions corresponding to conical surfaces of varying heights; this is valid even when the one-constant approximation is relaxed, see Sec. 1.3 in [44]. This gives us a subset of shapes described by $\zeta(r) = \pm m(r-1)$, with $m$ the vertical distance of the tip of the cone from the base, that allow us to move continuously from a flat disk to a cone of varying height. Minimizing Eq. (1) with respect to the height $m$ and for a fix phase $\psi$, one obtains the non-trivial condition

$$m^2 = \frac{A(1 - 2 \sin(\psi)^2) - 2 - B}{1 + A \sin(\psi)^2 + B}, \quad (2)$$

where $A = 4k/3k_B$ and $B = (4\sigma/k_B)(1 - \Delta^2)/\log(1/\Delta)$ are dimensionless parameters and $\Delta$ is a dimensionless cut-off lengthscale. The logarithmic divergence at $\Delta=0$ is commonly found in the context of topological defects [42] and also describes the divergence in curvature at the tip of the conical surface. We explore the effect of $\Delta$ on our results in Sec. 2 in [44]. Eq. (2) leads to the existence condition for conical shapes $A(1 - 2 \sin(\psi)^2) > 2 + B$. As observed, Eq. (2) shows that as $k/k_B$ or $k/\sigma$ increase, the minimal surface varies from a flat disc ($m=0$) to a cylinder ($m \to \infty$). The total free-energy associated with this minimal surface reads

$$\frac{\mathcal{F}_c}{\pi k_B} = \sqrt{(1 + A \sin(\psi)^2 + B)(A \cos(\psi)^2 - 1) \log(1/\Delta)}. \quad (3)$$

In addition, Eq. (3) shows that the minimal phase corresponds to $\psi=0$, even when the Frank constants $k_1=k_3$ are equal. The selection mechanism for the phase arises from the coupling between the director field and the extrinsic curvature, which tend to align the director field with the minimal principal curvature [18]. In the case of conical surfaces, the aster is favoured because it features only splay distortions, and vanishing twist and bend distortions. Furthermore, the total Frank free-energy for the aster decreases as the height of the cone increases, leading to the spontaneous out-of-plane deformation of a surface. This is balanced by the increased curvature and area of the surface. The threshold for a flat disc with an aster topological defect to become unstable is set when its energy (i.e. $2\pi k/3 \log(1/\Delta) + \sigma\pi(1 - \Delta^2)$) equals the energy (3) for $\psi=0$, that is when

$$k/3 = \frac{k_B}{2} + \sigma \frac{(1 - \Delta^2)}{\log(1/\Delta)}, \quad (4)$$

shown with the magenta line, Fig. 1b. This condition equals the existence condition of conical shapes for $\sigma=0$, Sec. 1.3 in [44]. A similar approach is used with varying elastic coefficients and defect phase to generate the green line in Fig. 1c.
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Height at the center and (b) phase at the boundary of the surface as a function of the three Frank constants with the constraint $k_1 + k_2 + k_3 = 1$. Magenta line is Eq. (4) and red line is Eq. (5). (c)&(d) Example surfaces with parameters: green star ($k_1 = 1/2$, $k_2 = 1/10$, $k_3 = 4/10$) blue star ($k_1 = 13/20$, $k_2 = 1/10$, $k_3 = 1/4$). (e) shows a conical aster surface while (d) features a concave surface with an aster at the center and spiral at the boundary. Blue (red) regions: negative (positive) Gaussian curvature. In all panels, $k_B = \sigma = 1/10$.</figcaption></figure>

According to Eq. 2, the height of the surface is also controlled by the phase of the defect. For equal Frank constants, the existence condition for conical shapes can be satisfied when $\psi < \pi/4$. Beyond this point, the increased energy associated with twist and bend distortions of the director field on a conical surface is too great and a flat surface is favoured. To further explore this effect, we numerically studied the equilibrium shapes of surfaces with a prescribed uniform phase, $\psi$, and a varying ratio of bend and splay elastic constants, $k_1/(k_1+k_3)$, Fig. 1c. In general, we find that vortices cannot deform a surface. In most cases, the aster generates the maximal out-of-plane deformation, except when $k_3 \gg k_1$ and the maximal deformation occurs near $\psi \sim \pi/8$, see Fig. 1e. This result is also predicted by the conical surface approximation, which is given by the green line in Fig. 1c.

To further explore these results, we now focus on the effects of varying all three Frank elastic constants; we also relax the constraint that $\psi$ is constant. Both $\zeta'$ and $\psi$ are only constrained by the rotational symmetry at the boundary. Figs. 2a and 2b show respectively the height at the center, $\zeta(0)$, and the phase at the outer boundary, $\psi(R)$, of the energy minimizing states for varying Frank coefficients. We found that the height of the deformed membrane is determined by the magnitude of $k_1$. Two transitions from flat to deformed surfaces were identified, which are primarily dependent on the relative values of $k_1$ and $k_3$. The transition from a flat to a deformed surface with an aster is determined by the threshold (4) with $k_1 = k/3$ (magenta line in Figs. 2a-b). A flat disc with a

<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) Height and (b) phase of the director field at the center of the membrane as a function of the phase at the boundary, $\psi(R) = \psi_R$, and the ratio between bend and splay coefficients ($k_2 = 1/3$) for surfaces with fix boundaries $\zeta'(R) = 0$. (c-e) Example of the three modes of out-of-plane deformation with fix boundary conditions: (c) pointy with an aster, (d) domed with an aster-to-vortex transition and (e) pointy with a vortex-to-aster transition. The states in (c-e) correspond to the black ($k_1 = 1/2$, $k_3 = 1/6$, $\psi_R = 0$), turquoise ($k_1 = 2/3$, $k_3 = 0$, $\psi_R = 0$) and pink ($k_1 = 0$, $k_3 = 2/3$, $\psi_R = \pi/2$) stars in (a) respectively. Blue (red) regions: negative (positive) Gaussian curvature. In all panels, $k_B = \sigma = 1/10$.</figcaption></figure>

vortex is linearly stable to out-of-plane deformations, see Sec. 1.4 in [44]. The red line in Fig. 2a-b can be found by comparing the energy of a flat disc with a vortex (i.e. $2\pi k_3 \log(1/\Delta) + \pi\sigma(1-\Delta^2)$) to the energy (3) for $\psi = 0$, that is when

$$ \frac{k_1}{k_B} = \frac{(k_3/k_B + B/4)^2}{(1+B)} + \frac{1}{4}. \qquad (5) $$

On a flat surface, if $k_1 > k_3$ ($k_1 < k_3$), the director field assumes a vortex (aster) configuration [42]. On a curved surface, however, an aster can be energetically favoured for values of $k_1 > k_3$; indeed all deformed surfaces here feature an aster at their core. In most cases, this is a conical aster deformation with constant phase, Fig. 2c. However, when $k_2 \ll k_3 < k_1$ we observe a new state with a spatially varying phase which features a conical aster close to the core combined with a negative Gaussian curvature spiral region near the boundary, Fig. 2d. The spiral region reduces splay at the cost of additional bend, however the bend is then exchanged for twist on the negative Gaussian curvature surface when $\psi$ has an intermediate value, which significantly reduces the free-energy density when $k_2$ is small.

Now, we consider the system under boundary conditions that might be found in biological systems. Hence, we fix $\psi(R) = \psi_R$ as the phase at the boundary of the membrane and fix $\zeta'(R) = 0$, implying that the surface must be flat at its boundary [45]. Figs. 3a&b show the
<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) Height at the apex as a function of the magnitude of $k_2$ for the three modes of out-of-plane deformation in Fig. 3c-e. The values of $k_1, k_3$ and $\psi_R$ are as outlined previously in caption Fig. 3. Here the constraint $k_1 + k_2 + k_3 = 1$ has been relaxed to allow $k_2$ to vary. (d-e) Examples of inverted versions for the three main deformation modes with $\zeta'(R) = 1$. Blue (red) regions: negative (positive) Gaussian curvature. In all panels, $k_B = \sigma = 1/10$.</figcaption></figure>

height and the phase at the center of the surface for a full range of the bend and splay constants ($k_2 = 1/3$) and the boundary phase $\psi_R$. On a flat surface, the value of $k_1/(k_1+k_3)$ controls the phase at the core of the topological defect. By varying $\psi_R$, we can explore regimes where the preferred phase at the center is frustrated with the boundary. We observe five distinct states. First, when the boundary phase $\psi_R$ is compatible with the dominant elastic coefficient, we observe flat asters and flat vortices. In addition to this there are three deformed configurations.

When $k_1 \approx k_3$ and $\psi_R \approx 0$ we find pointy surfaces featuring an aster, where the core and boundary of the defect are broadly in phase, see Fig. 3c. As in Fig. 2, the deformation here reduces splay distortions and hence stabilizes an aster, even in some cases when $k_1 > k_3$.

In the bend dominated regime $k_1 < k_3$ and for $\psi_R \approx \pi/2$, we observe pointy deformations on which the phase transitions from a vortex at the boundary to an aster at the center, Fig. 3e. In this configuration, the transition in the phase is localized to a ring of negative Gaussian curvature where bend distortions can be reduced at the cost of twist through the previously highlighted mechanism highlighted in Fig. 2d.

In the splay dominated regime $k_1 > k_3$ and for $\psi_R \approx 0$ we observe domed deformation on which the phase transitions from aster at the boundary to vortex at the center, Fig. 3d. Close to the boundary, the surface sharply deforms toward a cylindrical shape, abruptly reducing the splay at the cost of bend. On the cylindrical surface, the director field transitions from an aster to a vortex, further reducing splay and introducing additional twist and

bend distortions. The surface flattens toward the top featuring a vortex. Overall, this buckling mode is able to achieve the highest degree of deformation.

Twist distortions are only possible in spiral regions of the defect when $\psi$ has an intermediate value, this is only observed in the deformation modes that feature a transition in the phase. Therefore the value of $k_2$ can be used to mediate the deformation of these surfaces. The splay dominated deformation mode, Fig. 3d, features twist on an extended cylindrical section of the surface. When $k_2$ is increased the cylindrical area is reduced and the spiral is pushed into the positive Gaussian curvature, Fig. 4a&c. This reduces the induced twist and the height of the deformation reduces asymptotically toward a finite value, cyan curve Fig. 4a. The radial transition in the phase of the bend dominated deformation, Fig. 3e, is entirely within the negative Gaussian curvature region, which introduces increased twist. When $k_2$ is increased, the only way to reduce the twist is to reduce the negative Gaussian curvature, which in turn reduces the magnitude of the deformation and eventually suppress this mode, magenta curve Fig. 4a. Conversely, the conical aster deformation, Fig. 3c, features no transition in the phase, hence no twist, and is approximately unaffected by increasing $k_2$, black curve Fig. 4a. In fact, when $k_2$ is reduced to near zero, the director will locally distort to reduce splay and introduce twist, Fig. 4b.

We additionally observe bistability for these solutions. All surfaces shown thus far are identical under the reflection $\zeta \to -\zeta$. This symmetry can be broken by setting the boundary conditions $\zeta'(R) \neq 0$. All three previously identified buckling modes have a stable configuration in which $\zeta'$ changes sign closer to the center suggesting a similar mechanism could play a role in driving invagina- tion, see Figs. 4d-f.

Finally, this study reveals a novel mechanism for the spontaneous fusion of half-integer topological defects on fluid membranes. Consider the special case $\sigma = 0$ and $k_1 = k_2 = k_3 = k/3$. It is known that the energy minimizing nematic field with a total charge of +1 on a flat disc has two +1/2 defects. In this case, the energy scales, up to numerical pre-factors of order 1, as $F_f \sim \pi k/3 \log(1/\Delta)$ [46]. However, the free-energy of a conical surface with an aster at the apex scales sublinearly with $k$, (i.e. $F_c \sim \pi k_B \sqrt{(4k/3k_B - 1)} \log(1/\Delta)$). Therefore, a critical threshold $(k/3k_B)_c = \mathcal{O}(1)$ arises from the balance between the energies of these two states. Thus, if $(k/3k_B) > (k/3k_B)_c$, the pair of +1/2 topological defects can spontaneously fuse by deforming the surface out of plane. Indeed, past works have shown that this process can occur for prescribed conical surfaces [25].

This mechanism for spontaneous deformation of fluid membranes can drive out-of-plane deformations of biological system, such as a film of cytoskeletal filaments on a supported lipid bilayer or a cell monolayer on a supported elastic substrate, see Sec. 6 in [44]. For instance,
actin and microtubule films can exhibit nematic phases [47–49]. The measured reduced Frank constant is in the range of $k = 10^{-8} – 10^{-7} \mu N \mu m$ for actin filament films [50, 51] and is $k = 10^{-3} \mu N \mu m$ for microtubule films [52]. Using typical values for the bending rigidity and surface tension of lipid bilayers and ignoring other possible elastic contributions, allows one to evaluate the proximity to the instability threshold (4). In particular, in the bending dominate regime $k_B \gg \sigma R^2$, the critical threshold is of order of $(k/k_B)_c \sim 1$, and the ratio is $k/k_B = 1$ for a film of actin filaments and is $k/k_B = 10^4 – 10^5$ for a film of microtubules. In the tension dominated regime $k_B \ll \sigma R^2$, the critical threshold is $(k/\sigma R^2)_c \sim 1$ and therefore it can vary with the disc radius $R$. One can determine an upper bound by replacing the radius $R$ with the layer thickness $h$, which is the smallest lengthscale. This bound is $k/\sigma h^2 < 10$ for a film of both types of cytoskeletal filaments. These estimations suggest that in low-tension regimes, biological system can induce out-of-plane deformation by relieving elastic stresses at integer topological defects. In fact, this phenomenon may have been observed in Ref. [38] that studied a film of microtubules encapsulated on a vesicle. In this work, the authors reported that lowering the surface tension leads to spindle-like vesicles with two +1 defects localized at the spindle poles.

Shape dynamics driven by out-of-equilibrium processes could also be influenced by this mechanism. For instance, because the aster tends to facilitate surface deformations, it can be a nucleation point for out-of-plane deformations on active fluid membranes, like bacterial biofilms or cell monolayers [29, 37, 38, 40, 53].

We are grateful to Isabelle Bonnet, Mathieu Dedenon, Karsten Kruse, Jean François Joanny, Jacques Prost, and Feng-Ching Tsai for insightful discussions. DJGP acknowledges funding from the Swiss National Science Foundation under starting grant TMSGI2 211367. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 847718. Q.C acknowledges funding from Institut Curie EuReCa PhD Programme. This publication reflects only the author’s view and that the European Research Agency is not responsible for any use that may be made of the information it contains.

References cited in SI: [29, 37–40, 46–52, 54–76]

[1] A. Buguin, M.-H. Li, P. Silberzan, B. Ladoux, and P. Keller, Micro-actuators: When artificial muscles made of nematic liquid crystal elastomers meet soft lithography, Journal of the American Chemical Society **128**, 1088 (2006).

[2] T. H. Ware, M. E. McConney, J. J. Wie, V. P. Tondiglia, and T. J. White, Voxelated liquid crystal elastomers, Science **347**, 982 (2015).

[3] H. Aharoni, Y. Xia, X. Zhang, R. D. Kamien, and S. Yang, Universal inverse design of surfaces with thin nematic elastomer sheets, Proceedings of the National Academy of Sciences **115**, 7206 (2018).

[4] J. Pikul, S. Li, H. Bai, R. Hanlon, I. Cohen, and R. F. Shepherd, Stretchable surfaces with programmable 3d texture morphing for synthetic camouflaging skins, Science **358**, 210 (2017).

[5] E. Siéfert, E. Reyssat, J. Bico, and B. Roman, Bioinspired pneumatic shape-morphing elastomers, Nature materials **18**, 24 (2019).

[6] J. Dervaux and M. B. Amar, Morphogenesis of growing soft tissues, Physical review letters **101**, 068101 (2008).

[7] M. M. Müller, M. B. Amar, and J. Guven, Conical defects in growing sheets, Physical review letters **101**, 156104 (2008).

[8] C. D. Modes, K. Bhattacharya, and M. Warner, Disclination-mediated thermo-optical response in nematic glass sheets, Physical Review E **81**, 060701 (2010).

[9] C. D. Modes and M. Warner, Blueprinting nematic glass: Systematically constructing and combining active points of curvature for emergent morphology, Physical Review E **84**, 021711 (2011).

[10] E. Siéfert, E. Reyssat, J. Bico, and B. Roman, Programming stiff inflatable shells from planar patterned fabrics, Soft Matter **16**, 7898 (2020).

[11] C. Lee, M. Kim, Y. J. Kim, N. Hong, S. Ryu, H. J. Kim, and S. Kim, Soft robot review, International Journal of Control, Automation and Systems **15**, 3 (2017).

[12] O.-Y. Zhong-Can and W. Helfrich, Bending energy of vesicle membranes: General expressions for the first, second, and third variation of the shape energy and applications to spheres and cylinders, Physical Review A **39**, 5280 (1989).

[13] R. Lipowsky, The conformation of membranes, Nature **349**, 475 (1991).

[14] F. Jülicher and U. Seifert, Shape equations for axisymmetric vesicles: a clarification, Physical Review E **49**, 4728 (1994).

[15] J. Guven and P. Vázquez-Montejo, The geometry of fluid membranes: Variational principles, symmetries and conservation laws, The Role of Mechanics in the Study of Lipid Bilayers, 167 (2018).

[16] J. R. Frank and M. Kardar, Defects in nematic membranes can buckle into pseudospheres, Physical Review E **77**, 041705 (2008).

[17] G. Napoli and L. Vergori, Surface free energies for nematic shells, Physical Review E **85**, 061701 (2012).

[18] G. Napoli and L. Vergori, Extrinsic curvature effects on nematic shells, Physical review letters **108**, 207803 (2012).

[19] D. J. G. Pearce, Coupling the topological defect phase to the extrinsic curvature in nematic shells, Soft Matter **18**, 5082 (2022).

[20] A. Segatti, M. Snarski, and M. Veneroni, Equilibrium configurations of nematic liquid crystals on a torus, Physical Review E **90**, 012501 (2014).

[21] D. Jesenek, S. Kralj, R. Rosso, and E. G. Virga, Defect unbinding on a toroidal nematic shell, Soft matter **11**, 2434 (2015).

[22] I. Nitschke, M. Nestler, S. Praetorius, H. Löwen, and

<footnote>* daniel.pearce@unige.ch</footnote>

<footnote>† carles.blanch-mercader@curie.fr</footnote>
A. Voigt, Nematic liquid crystals on curved surfaces: a thin film limit, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences **474**, 20170686 (2018).

[23] M. Nestler, I. Nitschke, S. Praetorius, and A. Voigt, Orientational order on surfaces: The coupling of topology, geometry, and dynamics, Journal of Nonlinear Science **28**, 147 (2018).

[24] M. Nestler, I. Nitschke, H. Löwen, and A. Voigt, Properties of surface landau-de gennes q-tensor models, Soft Matter **16**, 4032 (2020).

[25] F. Vafa, G. H. Zhang, and D. R. Nelson, Defect absorption and emission for p-atic liquid crystals on cones, Physical Review E **106**, 024704 (2022).

[26] J. Santiago, Stresses in curved nematic membranes, Physical Review E **97**, 052706 (2018).

[27] J. Santiago, G. Chacón-Acosta, and F. Monroy, Membrane stress and torque induced by frank's nematic textures: A geometric perspective using surface-based constraints, Physical Review E **100**, 012704 (2019).

[28] D. J. G. Pearce, P. W. Ellis, A. Fernandez-Nieves, and L. Giomi, Geometrical control of active turbulence in curved topographies, Physical review letters **122**, 168002 (2019).

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

[30] D. Khoromskia and G. Salbreux, Active morphogenesis of patterned epithelial shells, eLife **12**, e75878 (2023).

[31] M. Nestler and A. Voigt, Active nematodynamics on curved surfaces—the influence of geometric forces on motion patterns of topological defects, arXiv preprint arXiv:2107.07779 (2021).

[32] L. A. Hoffmann, L. N. Carenza, J. Eckert, and L. Giomi, Theory of defect-mediated morphogenesis, Science Advances **8**, eabk2712 (2022).

[33] G. Salbreux, F. Jülicher, J. Prost, and A. Callan-Jones, Theory of nematic and polar active fluid surfaces, Physical Review Research **4**, 033158 (2022).

[34] S. Bell, S.-Z. Lin, J.-F. Rupprecht, and J. Prost, Active nematic flows over curved surfaces, Physical Review Letters **129**, 118001 (2022).

[35] F. Vafa, D. R. Nelson, and A. Doostmohammadi, Periodic orbits, pair nucleation, and unbinding of active nematic defects on cones, arXiv preprint arXiv:2310.06022 (2023).

[36] T. Singha, A. Polley, and M. Barma, Clustering of lipids driven by integrin, Soft Matter **19**, 6814 (2023).

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

[38] F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic, and A. R. Bausch, Topology and dynamics of active nematic vesicles, Science **345**, 1135 (2014).

[39] C. Blanch-Mercader, P. Guillamat, A. Roux, and K. Kruse, Quantifying material properties of cell monolayers by analyzing integer topological defects, Physical Review Letters **126**, 028101 (2021).

[40] Y. Ravichandran, M. Vogg, K. Kruse, D. J. Pearce, and A. Roux, Topology changes of the regenerating hydra define actin nematic defects as mechanical organizers of

morphogenesis, bioRxiv, 2024 (2024).

[41] F. C. Frank, I. liquid crystals. on the theory of liquid crystals, Discussions of the Faraday Society **25**, 19 (1958).

[42] P. De Gennes and J. Prost, The physics of liquid crystals, 1993, Oxford University Press, New York, Olbrich E., Marinov O., Davidov D., Phys. Rev. E **2713**, 48 (1993).

[43] M. Nestler, I. Nitschke, and A. Voigt, A finite element approach for vector-and tensor-valued surface pdes, Journal of Computational Physics **389**, 48 (2019).

[44] Supplementary material, URL **will_be_inserted_byPublisher**.

[45] Note that the boundary condition $\zeta'(R)=0$ enforces that the total Gaussian curvature, including the tip, is zero.

[46] G. Duclos, C. Erlenkämper, J.-F. Joanny, and P. Silberzan, Topological defects in confined populations of spindle-shaped cells, Nature Physics **13**, 58 (2017).

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

[48] A. Sciortino and A. R. Bausch, Pattern formation and polarity sorting of driven actin filaments on lipid membranes, Proceedings of the National Academy of Sciences **118**, e2017047118 (2021).

[49] F. L. Memarian, J. D. Lopes, F. J. Schwarzendahl, M. G. Athani, N. Sarpangala, A. Gopinathan, D. A. Beller, K. Dasbiswas, and L. S. Hirst, Active nematic order and dynamic lane formation of microtubules driven by membrane-bound diffusing motors, Proceedings of the National Academy of Sciences **118**, e2117107118 (2021).

[50] R. Zhang, N. Kumar, J. L. Ross, M. L. Gardel, and J. J. De Pablo, Interplay of structure, elasticity, and dynamics in actin-based nematic materials, Proceedings of the National Academy of Sciences **115**, E124 (2018).

[51] V. Yadav, D. S. Banerjee, A. P. Tabatabai, D. R. Kozvar, T. Kim, S. Banerjee, and M. P. Murrell, Filament nucleation tunes mechanical memory in active polymer networks, Advanced functional materials **29**, 1905243 (2019).

[52] I. Vélez-Cerón, P. Guillamat, F. Sagués, and J. Ignés-Mullol, Probing active nematics with in-situ microfabricated elastic inclusions, arXiv preprint arXiv:2307.11587 (2023).

[53] M. Prasad, N. Obana, S.-Z. Lin, S. Zhao, K. Sakai, C. Blanch-Mercader, J. Prost, N. Nomura, J.-F. Rupprecht, J. Fattaccioli, et al., Alcanivorax borkumensis biofilms enhance oil degradation by interfacial tubulation, Science **381**, 748 (2023).

[54] K. D. Endresen, M. Kim, M. Pittman, Y. Chen, and F. Serra, Topological defects of integer charge in cell monolayers, Soft Matter **17**, 5878 (2021).

[55] M. J. Bowick and L. Giomi, Two-dimensional matter: order, curvature and defects, Advances in Physics **58**, 449 (2009).

[56] E. Evans and D. Needham, Physical properties of surfactant bilayer membranes: thermal transitions, elasticity, rigidity, cohesion and colloidal interactions, Journal of Physical Chemistry **91**, 4219 (1987).

[57] W. Rawicz, K. C. Olbrich, T. McIntosh, D. Needham, and E. Evans, Effect of chain length and unsaturation on elasticity of lipid bilayers, Biophysical journal **79**, 328 (2000).

[58] A. Roux, G. Koster, M. Lenz, B. Sorre, J.-B. Manneville, P. Nassoy, and P. Bassereau, Membrane curvature controls dynamin polymerization, Proceedings of the Na-
tional Academy of Sciences **107**, 4141 (2010).

[59] R. Phillips, T. Ursell, P. Wiggins, and P. Sens, Emerging roles for lipids in shaping membrane-protein function, Nature **459**, 379 (2009).

[60] F. Hanakam, R. Albrecht, C. Eckerskorn, M. Matzner, and G. Gerisch, Myristoylated and non-myristoylated forms of the ph sensor protein hisactophilin ii: intracellular shuttling to plasma membrane and nucleus monitored in real time by a fusion with green fluorescent protein., The EMBO journal **15**, 2935 (1996).

[61] A. G. Clark, K. Dierkes, and E. K. Paluch, Monitoring actin cortex thickness in live cells, Biophysical journal **105**, 570 (2013).

[62] V. Laplaud, N. Levernier, J. Pineau, M. S. Roman, L. Barbier, P. J. Sáez, A.-M. Lennon-Duménil, P. Vargas, K. Kruse, O. Du Roure, et al., Pinching the cortex of live cells reveals thickness instabilities caused by myosin ii motors, Science Advances **7**, eabe3640 (2021).

[63] D. Oriola, F. Jülicher, and J. Brugués, Active forces shape the metaphase spindle through a mechanical instability, Proceedings of the National Academy of Sciences **117**, 16154 (2020).

[64] A. R. Harris, L. Peter, J. Bellis, B. Baum, A. J. Kabla, and G. T. Charras, Characterizing the mechanics of cultured cell monolayers, Proceedings of the National Academy of Sciences **109**, 16449 (2012).

[65] A. D. S. Da Costa, R. Subbiah, S. J. Oh, H. Jeong, J.-I. Na, K. Park, I.-S. Choi, and J. H. Shin, Fibroblasts close a void in free space by a purse-string mechanism, ACS applied materials & interfaces **14**, 40522 (2022).

[66] E. Latorre, S. Kale, L. Casares, M. Gómez-González, M. Uroz, L. Valon, R. V. Nair, E. Garreta, N. Montserrat, A. Del Campo, et al., Active superelasticity in three-dimensional epithelia of controlled shape, Nature **563**, 203 (2018).

[67] A. Trushko, I. Di Meglio, A. Merzouki, C. Blanch-Mercader, S. Abuhattum, J. Guck, K. Alessandri, P. Nasso, K. Kruse, B. Chopard, et al., Buckling of an epithelium growing under spherical confinement, Developmental cell **54**, 655 (2020).

[68] J. Fouchard, T. P. Wyatt, A. Proag, A. Lisica, N. Khalilgharibi, P. Recho, M. Suzanne, A. Kabla, and G. Charras, Curling of epithelial monolayers reveals coupling between active bending and tissue tension, Proceedings of the National Academy of Sciences **117**, 9377 (2020).

[69] J. Duque, A. Bonfanti, J. Fouchard, E. Ferber, A. Harris, A. Kabla, and G. Charras, Fracture in living cell monolayers, BioRxiv, 2023 (2023).

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

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

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

[73] P. Lee and C. W. Wolgemuth, Crawling cells can close wounds without purse strings or signaling, PLoS computational biology **7**, e1002007 (2011).

[74] C. Pérez-González, R. Alert, C. Blanch-Mercader, M. Gómez-González, T. Kolodziej, E. Bazellieres, J. Casademunt, and X. Trepat, Active wetting of epithe-

lial tissues, Nature physics **15**, 79 (2019).

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

[76] P. Guillamat, J. Ignés-Mullol, S. Shankar, M. C. Marchetti, and F. Sagués, Probing the shear viscosity of an active nematic film, Physical review E **94**, 060602 (2016).
## SUPPLEMENTAL MATERIAL:

# SECTION 1: THEORETICAL DESCRIPTION

We consider a surface with the following free-energy

$$ \mathcal{F} = \int_{\mathcal{A}} \{k_B \mathcal{H}^2 + k_1 (\nabla \cdot \hat{n})^2 + k_2 (\hat{n} \cdot (\nabla \times \hat{n}))^2 + k_3 (\hat{n} \times (\nabla \times \hat{n}))^2 + \sigma\} da. \quad (S1) $$

The first term is the bending energy with mean curvature $\mathcal{H}$, the second, third and fourth terms are the Frank free-energy associated respectively with splay, twist and bend distortions of the director field $\hat{n}$ and the fifth term is the compressional energy. In the same order, the corresponding elastic coefficients are: the bending rigidity $k_B$, the reduced Frank constant $k_1$, $k_2$, and $k_3$ and the surface tension $\sigma$. In the following, we assume that the system is deep into the nematic phase and impose that $|\hat{n}| = 1$.

### Section 1.1: Lagrangian equations for surfaces of revolution with an integer topological defect

In this subsection, we derive the Lagrangian equations for a surface of revolution with an integer topological defect at the axis of rotation.

Being $\mathbf{r}$ the position vector, the surface is parametrised as follows,

$$ \mathbf{r} = (r(s) \cos(\theta), r(s) \sin(\theta), \zeta(s)) \quad (S2) $$

where $\theta$ is the azimuthal angle and $s$ is the arc-length and it satisfies the constrain $|\partial_s \mathbf{r}| = 1$. Because of the previous constrain, one can define $\partial_s r = \cos(\phi(s))$ and $\partial_s \zeta = -\sin(\phi(s))$, where $\phi(s)$ is minus the angle of the tangent vector of the generatrix with respect to the radial direction. In addition, because of the symmetries of surfaces of revolution, we ignore variations in the $\theta$ coordinate. Hence, from now on, we denote by primes the derivatives with respect to the arclength $s$. The variation of the position vector (S2) with respect to $\theta$ and $s$ defines an orthogonal base on the tangent plane given by

$$ \partial_s \mathbf{r} = (r'(s) \cos(\theta), r'(s) \sin(\theta), \zeta'(s)) \quad (S3a) $$

$$ \partial_{\theta} \mathbf{r} = r(s)(-\sin(\theta), \cos(\theta), 0) \quad (S3b) $$

which can be normalised as $\hat{e}_s = \partial_s \mathbf{r}$ and $\hat{e}_{\theta} = \partial_{\theta} \mathbf{r}/r(s)$. A normal to the tangent plane reads

$$ \hat{n} = (-\zeta'(s) \cos(\theta), -\zeta'(s) \sin(\theta), r'(s)) \quad (S4) $$

The components of the metric tensor $g_{\alpha\beta}$ read

$$ g_{ss} = \partial_s \mathbf{r} \cdot \partial_s \mathbf{r} = 1 \quad g_{\theta\theta} = r^2 \quad g_{s\theta} = g_{\theta s} = 0 \quad (S5) $$

and the differential of area on the surface reads $da = \sqrt{\det(g_{\alpha\beta})} ds d\theta = r ds d\theta$.

The components of the curvature tensor read

$$ C_{ss} = -\partial_s^2 \mathbf{r} \cdot \hat{n} = \phi' \quad C_{\theta\theta} = -\partial_{\theta\theta}^2 \mathbf{r} \cdot \hat{n} = -\zeta' r \quad C_{s\theta} = C_{\theta s} = -\partial_{s\theta}^2 \mathbf{r} \cdot \hat{n} = 0 \quad (S6) $$

and the mean curvature as $\mathcal{H} = (C_{ss}/g_{ss} + C_{\theta\theta}/g_{\theta\theta})/2 = (\phi' - \zeta'/r)/2$.

The variation of our tangent plane basis with respect to the parameters $\theta$ and $s$ reads

$$ \partial_s \hat{e}_s = -\phi' \hat{n} \quad \partial_{\theta} \hat{e}_s = r' \hat{e}_{\theta} \quad \partial_s \hat{e}_{\theta} = 0 \quad \partial_{\theta} \hat{e}_{\theta} = \zeta' \hat{n} - r' \hat{e}_s \quad (S7) $$

which allow one to compute the tensor of distortions of the director field $\partial_{\alpha} \hat{n}_{\beta}$. We assume that the director field is prescribed and corresponds to that of an integer defect at the axis of rotation. In addition, we assume that the director field is constrained to be parallel to the surface of revolution (i.e. $\hat{n} \cdot \hat{n} = 0$). Therefore $\hat{n} = \cos(\psi) \hat{e}_s + \sin(\psi) \hat{e}_{\theta}$, where $\psi$ is the angle of the director field with respect to the direction $\hat{e}_s$. For simplicity, $\psi$ is considered uniform. In this
case, $\psi = 0$ corresponds to an aster, $0 < \psi < \pi/2$ to a spiral, and $\psi = \pi/2$ to a vortex. Under these approximations and using the relations (S7), the divergence and the curl of the director field reads

$$ \nabla \cdot \hat{n} = (\hat{e}_s \partial_s + \frac{\hat{e}_\theta}{r} \partial_\theta) \cdot (\cos(\psi) \hat{e}_s + \sin(\psi) \hat{e}_\theta) \quad (S8a) $$

$$ = \cos(\psi)\hat{e}_s \cdot \partial_s \hat{e}_s + \sin(\psi)\hat{e}_s \cdot \partial_s \hat{e}_\theta + \frac{\cos(\psi)}{r}\hat{e}_\theta \cdot \partial_\theta \hat{e}_s + \frac{\sin(\psi)}{r}\hat{e}_\theta \cdot \partial_\theta \hat{e}_\theta \quad (S8b) $$

$$ = \frac{\cos(\psi)r'}{r} \quad (S8c) $$

$$ \nabla \times \hat{n} = (\hat{e}_s \partial_s + \frac{\hat{e}_\theta}{r} \partial_\theta) \times (\cos(\psi) \hat{e}_s + \sin(\psi) \hat{e}_\theta) \quad (S8d) $$

$$ = \cos(\psi)\hat{e}_s \times \partial_s \hat{e}_s + \sin(\psi)\hat{e}_s \times \partial_s \hat{e}_\theta + \frac{\cos(\psi)}{r}\hat{e}_\theta \times \partial_\theta \hat{e}_s + \frac{\sin(\psi)}{r}\hat{e}_\theta \times \partial_\theta \hat{e}_\theta \quad (S8e) $$

$$ = \phi' \cos(\psi) \hat{e}_\theta + \frac{\sin(\psi)\zeta'}{r} \hat{e}_s + \frac{\sin(\psi)r'}{r} \hat{n} \quad (S8f) $$

Therefore the splay, twist and bend distortions of the director field in the free-energy (S1) read,

$$ k_1(\nabla \cdot \hat{n})^2 = k_1 \left( \frac{\cos(\psi)r'}{r} \right)^2 \quad (S9a) $$

$$ k_2(\hat{n} \cdot (\nabla \times \hat{n}))^2 = k_2 \left( \left( \phi' + \frac{\zeta'}{r} \right) \cos(\psi) \sin(\psi) \right)^2 \quad (S9b) $$

$$ k_3(\hat{n} \times (\nabla \times \hat{n}))^2 = k_3 \left( \left( \frac{\sin(\psi)r'}{r} \right)^2 + \left( \phi' \cos(\psi)^2 - \frac{\sin(\psi)^2 \zeta'}{r} \right)^2 \right) \quad (S9c) $$

In this case, the free-energy (S1) takes the form

$$ \begin{aligned} \mathcal{F} = & \int_{\mathcal{A}} \left\{ \frac{k_B}{4} \left( \phi' - \frac{\zeta'}{r} \right)^2 + k_1 \left( \frac{\cos(\psi)r'}{r} \right)^2 + k_2 \left( \left( \phi' + \frac{\zeta'}{r} \right) \cos(\psi) \sin(\psi) \right)^2 \right. \\ & \left. + k_3 \left( \left( \frac{\sin(\psi)r'}{r} \right)^2 + \left( \phi' \cos(\psi)^2 - \frac{\sin(\psi)^2 \zeta'}{r} \right)^2 \right) + \sigma \right\} r d\theta ds \\ & + \int_{\mathcal{A}} \gamma(r' - \cos(\phi)) + \eta(\zeta' + \sin(\phi)) ds d\theta = 2\pi \int_{\mathcal{A}} \mathcal{L}[\phi, \phi', \zeta', r, r'] ds \end{aligned} \quad (S10) $$

where $\gamma$ and $\eta$ are two Lagrange multipliers to enforce the constraints associated with the arclength parametrisation.

From Eq. (S10), one can identify a Lagrangian $\mathcal{L}$. Then the Lagrangian equations for variations of $\phi$, $r$ and $\zeta$ are given by

$$ \delta\phi : \frac{\partial\mathcal{L}}{\partial\phi} - \frac{d}{ds}\left(\frac{\partial\mathcal{L}}{\partial\phi'}\right) = 0 \quad (S11a) $$

$$ \delta r : \frac{\partial \mathcal{L}}{\partial r} - \frac{d}{ds} \left( \frac{\partial \mathcal{L}}{\partial r'} \right) = 0 \quad (S11b) $$

$$ \delta \zeta : \frac{\partial \mathcal{L}}{\partial \zeta} - \frac{d}{ds} \left( \frac{\partial \mathcal{L}}{\partial \zeta'} \right) = 0 \quad (S11c) $$

with the boundary contributions

$$ \delta\phi : \frac{\partial\mathcal{L}}{\partial\phi'} = 0 \quad (S12a) $$

$$ \delta r : \frac{\partial \mathcal{L}}{\partial r'} = 0 \quad (S12b) $$

$$ \delta \zeta : \frac{\partial \mathcal{L}}{\partial \zeta'} = 0 \quad (S12c) $$

Because the Lagrangian (S10) is independent on the arclength parameter $s$, one can identify a Hamiltonian $H$, which is equal to a constant. When the total length is allowed to vary during the energy minimisation process, then this constant is zero. We provide the explicit forms of Eqs. (S11), (S12) and the Hamiltonian in Appendix 1.
## Section 1.2: Normal force balance

In this section, we compute the normal force balance equation for a surface of revolution with a prescribed integer topological defect as described in Section 1.1.

To obtain the normal force balance equation, we combine Eqs. (S30) and (S32). First, using Eq. (S30b) and Eq. (S32), one can obtain

$$ (\bar{\eta}r' - \gamma\zeta')' = -\gamma'\zeta' + \phi'(\bar{\eta}\zeta' + \gamma r') \\ = -\left\{\frac{k_B}{4}\left(\phi'^2 - \left(\frac{\zeta'}{r}\right)^2\right) + (\phi' \cos(\psi))^2(k_2 \sin(\psi)^2 + k_3 \cos(\psi)^2)\right. \\ \left. - \left(\frac{\zeta' \sin(\psi)}{r}\right)^2(k_2 \cos(\psi)^2 + k_3 \sin(\psi)^2) - (k_1 \cos(\psi)^2 + k_3 \sin(\psi)^2)\left(\frac{r'}{r}\right)^2\right\} \left(\phi' + \frac{\zeta'}{r}\right)r \\ + \sigma r \left(\phi' - \frac{\zeta'}{r}\right) + 2(k_1 \cos(\psi)^2 + k_3 \sin(\psi)^2)\left(\frac{r'\zeta'}{r}\right)' \quad (S13a) $$

where we used the condition $\bar{\eta}' = 0$ given by Eq. (S30c). This expression allows one to re-express the derivative with respect to the arclength of Eq. (S30a) as

$$ 0 = -\left\{\frac{k_B}{4}\left(\phi'^2 - \left(\frac{\zeta'}{r}\right)^2\right) + (\phi' \cos(\psi))^2(k_2 \sin(\psi)^2 + k_3 \cos(\psi)^2)\right. \\ \left. - \left(\frac{\zeta' \sin(\psi)}{r}\right)^2(k_2 \cos(\psi)^2 + k_3 \sin(\psi)^2) - (k_1 \cos(\psi)^2 + k_3 \sin(\psi)^2)\left(\frac{r'}{r}\right)^2\right\} \left(\phi' + \frac{\zeta'}{r}\right)r \\ + \sigma r \left(\phi' - \frac{\zeta'}{r}\right) - \left(\frac{r'\zeta'}{r}\right)' \left(\frac{k_B}{2} + 2 \sin(\psi)^2(k_2 \cos(\psi)^2 + k_3 \sin(\psi)^2) - 2(k_1 \cos(\psi)^2 + k_3 \sin(\psi)^2)\right) \\ \left. - (\phi'r)'' \left(\frac{k_B}{2} + 2 \cos(\psi)^2(k_2 \sin(\psi)^2 + k_3 \cos(\psi)^2)\right)\right\} \quad (S14) $$

Up to a scaling factor $r$, Eq. (S14) is the normal force balance. Indeed, the terms proportional to the bending rigidity $k_B$ and the terms proportional to the surface tension $\sigma$ are equivalent to the normal forces found by O-Y Zhong-can and W. Helfrich. Similarly, the terms proportional to the splay, twist and bend Frank elastic constant $k_1$, $k_2$, and $k_3$ are equivalent to the normal force found by J.A. Santiago et al (2019).

## Section 1.3: Exact solutions of the normal force balance equation

In this section, we compute exact solution of the normal force balance equation and analyse their stability conditions. In the special case where $\phi = \phi_0$ is constant and $\sigma = 0$, equation (S14) reduces to

$$ 0 = -r \left(\frac{\sin(\phi_0)}{r}\right) \left\{ \left(\frac{\sin(\phi_0)}{r}\right)^2 \left(\frac{k_B}{4} + \sin(\psi)^2(k_2 \cos(\psi)^2 + k_3 \sin(\psi)^2)\right) + \left(\frac{\cos(\phi_0)}{r}\right)^2 \left(\frac{k_B}{2} + 2 \sin(\psi)^2(k_2 \cos(\psi)^2 + k_3 \sin(\psi)^2) - (k_1 \cos(\psi)^2 + k_3 \sin(\psi)^2)\right) \right\} \quad (S15) $$

which admits three solutions: the trivial solution $\phi_0 = 0$ corresponding to a flat disc, and two non-trivial solutions corresponding to cones with an inclination angle (complementary to the opening angle) given by

$$ \tan(\phi_0)^2 = \frac{(k_1 \cos(\psi)^2 + k_3 \sin(\psi)^2) - \frac{k_B}{2} - 2 \sin(\psi)^2(k_2 \cos(\psi)^2 + k_3 \sin(\psi)^2)}{\frac{k_B}{4} + \sin(\psi)^2(k_2 \cos(\psi)^2 + k_3 \sin(\psi)^2)} \quad (S16) $$

The existence condition for these non-trivial solutions is that

$$ (k_1 \cos(\psi)^2 + k_3 \sin(\psi)^2) > \frac{k_B}{2} + 2 \sin(\psi)^2(k_2 \cos(\psi)^2 + k_3 \sin(\psi)^2) \quad (S17) $$
For an aster $\psi = 0$, steady-state cones exists when $k_1 > k_B/2$, whereas for a vortex $\psi = \pi/2$, no steady-state cone are expected. The existence of steady-state conical surfaces, for intermediate spirals $0 < \psi < \pi/2$, depends on the values of the elastic constants.

Next we analyse the stability of these steady-state solutions by computing their total free-energy. Considering that $\phi = \phi_0$ is constant and $\sigma = 0$, Eq. (S10) reduces to

$$ \mathcal{F} = 2\pi \left\{ (\zeta')^2 \left( \frac{k_B}{4} + \sin(\psi)^2(k_2 \cos(\psi)^2 + k_3 \sin(\psi)^2) \right) + (r')^2 \left( k_1 \cos(\psi)^2 + k_3 \sin(\psi)^2 \right) \right\} \frac{\log(R/a)}{r'} \quad (S18) $$

where we used that $dr = r'ds$, $\zeta' = -\sin(\phi_0)$ and $r' = \cos(\phi_0)$, and considered as boundary conditions that $r(s_2) = R$ and $r(s_1) = a$, where R is the radius of the disc and a is a microscopic length. Eq. (S18) corresponds to the total free-energy of a integer defect with a constant angle $\psi$ at the tip of a conical surface with an inclination angle $\phi_0$.

Minimization of Eq. (S18) with respect to $\phi_0$ leads to same three solutions found by solving the force balance equation. Minimization of Eq. (S18) with respect to $\psi$ leads to three types of topological defects: an aster $\psi = 0$, a vortex $\psi = \pi/2$, and a spiral with an angle that satisfies

$$ \tan(\phi_0)^2 = \frac{k_1 - k_3}{2k_3 \sin(\psi)^2 + k_2(\cos(\psi)^2 - \sin(\psi)^2)} \quad (S19) $$

The minimal surfaces with an aster defect are the flat disc with total free-energy

$$ \mathcal{F}/2\pi = k_1 \log(R/a) \quad (S20) $$

and the conical surface with elevation angle and total free-energy

$$ \tan(\phi_0)^2 = \frac{k_1 - k_B/2}{k_B/4} > 0 \quad (S21) $$

$$ \mathcal{F}/2\pi = 2\sqrt{(k_1 - k_B/4)k_B/4 \log(R/a)} \quad (S22) $$

The minimal surface with a vortex defect is a flat disc with total free-energy

$$ \mathcal{F}/2\pi = k_3 \log(R/a) \quad (S23) $$

Finally, the minimal surface with the spiral defect that satisfies (S19) is a conical surface (except for $k_1 = k_3$) with the elevation angle (??) and total free-energy (S18).

Stability conditions can be identified by comparing the energies of an aster on a minimal conical surface, an aster on a flat disc and a vortex on a flat disc. Respectively, these conditions read

$$ \frac{k_1}{k_B} = \frac{1}{2} \quad (S24) $$

$$ \frac{k_1}{k_B} = \left( \frac{k_3}{k_B} \right)^2 + \frac{1}{4} \quad \text{with} \quad \frac{k_1}{k_B} > \frac{1}{2} \quad (S25) $$

$$ \frac{k_1}{k_B} = \frac{k_3}{k_B} \quad (S26) $$

## Section 1.4: Linear stability of a flat disc with an integer topological defect to shape fluctuations

In this section, we perform a linear stability analysis of a flat disc with an integer topological defect. Next, we discuss the parametric conditions whereby this state can become unstable as a function of the type of topological defect.

As shown above, the flat configuration is a steady state solution of the normal force balance equation (S14). Here, we study its stability to small amplitude shape fluctuations. Therefore we consider $\phi = \delta\phi \ll 1$ and its derivatives to be small. Note that $r' = 1 + \mathcal{O}(\delta\phi^2)$, and $\zeta' = -\delta\phi + \mathcal{O}(\delta\phi^2)$.
To linear order in perturbations $\delta\phi$, the Lagrangian equations (S30b) and (S30c) provide a solution for the Lagrange multipliers $\gamma$ and $\bar{\eta}$. We enforce that their integration constant vanish. Inserting their expression into (S30a), one obtains the equation

$$0 = A(r) \frac{\phi}{r} - B(\phi'r)' \qquad (S27)$$

where $A = -(k_1 \cos(\psi)^2 + k_3 \sin(\psi)^2) + \sigma r^2 + k_B/2 + 2 \sin(\psi)^2(k_2 \cos(\psi)^2 + k_3 \sin(\psi)^2)$ and $B = k_B/2 + 2 \cos(\psi)^2(k_2 \sin(\psi)^2 + k_3 \cos(\psi)^2) > 0$ is a positive coefficient.

In the special case $\sigma = 0$, the general solution of Eq. (S27) reads

$$\phi = C \cosh \left( \sqrt{A/B} \log(r) + D \right) \qquad (S28)$$

where $C$ and $D$ are integration constants to be determined. Note that the solution (S28) has a logarithmic divergence at $r \to 0$ that can not be eliminated by a specific choice of the integration constants. When $A < 0$, the solution (S28) becomes oscillatory which is a signature of an instability. Therefore, we identify the threshold whereby a flat disc with an integer topological defect can become unstable as $A = 0$, or equivalently

$$k_1 \cos(\psi)^2 + k_3(\sin(\psi)^2 - 2 \sin(\psi)^4) = 2k_2 \sin(\psi)^2 \cos(\psi)^2 + k_B/2 \qquad (S29)$$

Note that this condition matches with the existence condition of conical surfaces (S17) found for $\sigma = 0$. Equation (S29) reveals two linear mechanisms for generating out-of-plane deformations. On the one hand, in the splay dominated regime $k_1 \gg k_2, k_3$, the first term on the left hand side of the stability condition (S29) dominates. The critical splay Frank constant at which the flat disc becomes linearly unstable increase as the phase of the defect $\psi$ increases, and diverges for a vortex with $\psi = \pi/2$. On the other hand, in the bend dominated regime $k_3 \gg k_1, k_2$, the stability condition (S29) can be satisfied for a range of phase $\psi$ from 0 to $\pi/4$. Therefore, topological defects with a phase $\psi > \pi/4$ are unable to linearly destabilize a flat disc by decreasing bend distortions. Finally, in the twist dominated regime $k_2 \gg k_1, k_3$, flat discs are linearly stable irrespective of the type of topological defect. Therefore, we conclude that the aster and all spirals, except for the vortex, are able to generate spontaneous out-of-plane deformations on a flat disc via two linear mechanisms that reduce the energy associated with either splay or bend distortions by deforming the surface in the out-of-plane direction.

In the general case $\sigma \neq 0$, we estimated an approximate instability threshold by comparing the total energies of a conical surface with an aster to the total energy of a flat disc with either an aster or a vortex. In this case, we find that the instability threshold is increased by the stabilization of flat geometries by surface tension.

## Appendix 1: Lagrangian equations and Hamiltonian

Then the equilibrium equations for variations of $\phi$, $r$ and $\zeta$ are given by

$$\begin{aligned} \delta\phi: 0 = \bar{\eta}r' - \gamma\zeta' - \frac{\zeta'r'}{r} \left( \frac{k_B}{2} + 2\sin(\psi)^2(k_2\cos(\psi)^2 + k_3\sin(\psi)^2) \right) \\ - (\phi'r)' \left( \frac{k_B}{2} + 2\cos(\psi)^2(k_2\sin(\psi)^2 + k_3\cos(\psi)^2) \right) \end{aligned} \qquad (S30a)$$

$$\begin{aligned} \delta r: 0 = \phi'^2 \left( \frac{k_B}{4} + \cos(\psi)^2(k_2 \sin(\psi)^2 + k_3 \cos(\psi)^2) \right) + \left( \left( \frac{r'}{r} \right)^2 - \frac{2r''}{r} \right) (k_1 \cos(\psi)^2 + k_3 \sin(\psi)^2) \\ - \gamma' + \sigma - \left( \frac{\zeta'}{r} \right)^2 \left( \frac{k_B}{4} + \sin(\psi)^2(k_2 \cos(\psi)^2 + k_3 \sin(\psi)^2) \right) \end{aligned} \qquad (S30b)$$

$$\delta\zeta: 0 = (\bar{\eta})' \qquad (S30c)$$

with the boundary conditions

$$\delta\phi: 0 = \phi'r \left( \frac{k_B}{2} + 2\cos(\psi)^2(k_2\sin(\psi)^2 + k_3\cos(\psi)^2) \right) - \zeta' \left( \frac{k_B}{2} + 2(k_3 - k_2)\cos(\psi)^2\sin(\psi)^2 \right) \qquad (S31a)$$

$$\delta r: 0 = \frac{2r'}{r} \left( k_1 \cos(\psi)^2 + k_3 \sin(\psi)^2 \right) + \gamma \qquad (S31b)$$

$$\delta\zeta: 0 = \bar{\eta} \qquad (S31c)$$
where we defined $\bar{\eta} = \eta - \phi'(k_B/2 + 2(k_3 - k_2)\sin(\psi)^2\cos(\psi)^2) + \zeta'(k_B/2 + 2\sin(\psi)^2(k_2\cos(\psi)^2 + k_3\sin(\psi)^2))/r$.
The Hamiltonian reads

$$ H = -\mathcal{L} + \phi' \frac{\partial \mathcal{L}}{\partial \phi'} + r' \frac{\partial \mathcal{L}}{\partial r'} + \zeta' \frac{\partial \mathcal{L}}{\partial \zeta'} \\ = \frac{rk_B}{4} \left( \phi'^2 - \left(\frac{\zeta'}{r}\right)^2 \right) + rk_1 \left( \frac{r' \cos(\psi)}{r} \right)^2 + rk_2 \cos(\psi)^2 \sin(\psi)^2 \left( \phi'^2 - \left(\frac{\zeta'}{r}\right)^2 \right) \\ + rk_3 \left( \left( \frac{r' \sin(\psi)}{r} \right)^2 + (\phi' \cos(\psi))^2 - \left( \frac{\zeta' \sin(\psi)}{r} \right)^2 \right) - \sigma r + \bar{\eta} \zeta' + \gamma r' \quad (S32) $$

# SECTION 2: HEIGHT OF A CONICAL SURFACE WITH ARBITRARY PHASE

## Section 2.1: Energy of a conical surface

We assume a cylindrical Monge gauge with $\zeta = mr$ and $\psi$ constant. The mean curvature of such a surface is given by $\mathcal{H} = \frac{m}{2r\sqrt{1+m^2}}$, the area element of the surface is given by $r\sqrt{1+m^2}$. We take the radius of the membrane to be $R=1$.

With this we can write down the free-energy associated with the shape of the surface as

$$ F_s = 2\pi \int_{\Delta}^{1} [k_B \mathcal{H}^2 + \sigma] r \sqrt{1+m^2} dr \quad (S33) $$

$$ = \frac{\pi k_B m^2}{2\sqrt{1+m^2}} \log(1/\Delta) + \sigma\pi\sqrt{1+m^2}[1-\Delta^2], \quad (S34) $$

where we have introduced $\Delta$ as the inner radius of the surface. This is to account for the fact that on a cone $\mathcal{H}$ diverges at $r=0$. For comparison with numerical results, $\Delta = 1/N$, where $N$ is the number of discretization bins used to evaluate the functions.

The energy associated with distortions in the nematic director is given by

$$ F_{LC} = 2\pi \int_{\Delta}^{1} \frac{k_1 \cos^2(\psi)}{r\sqrt{1+m^2}} + \frac{k_2 m^2 \sin^2(\psi) \cos^2(\psi)}{r\sqrt{1+m^2}} + \frac{k_3 \sin^2(\psi)(1+m^2 \sin^2(\psi))}{r\sqrt{1+m^2}} dr \quad (S35) $$

$$ = \frac{2\pi \log(1/\Delta)}{\sqrt{1+m^2}} \left[ k_1 \cos^2(\psi) + k_2 m^2 \sin^2(\psi) \cos^2(\psi) + k_3 \sin^2(\psi)(1+m^2 \sin^2(\psi)) \right] \quad (S36) $$

We introduce $L = 2\pi \log(1/\Delta)$ and $M = \pi[1 - \Delta^2]$. We now write the total free-energy of a conical surface with a defect of constant phase $\psi$ as

$$ F = L \left[ k_1 \cos^2(\psi) + k_3 \sin^2(\psi) \right] / \sqrt{1+m^2} + L \left[ k_2 \sin^2(\psi) \cos^2(\psi) + k_3 \sin^4(\psi) + k_B/4 \right] m^2 / \sqrt{1+m^2} + M\sigma \sqrt{1+m^2} \quad (S37) $$

## Section 2.2: Equilibrium height of a conical surface

Taking the derivative with respect to $m$ and setting it to zero, we arrive at an expression for the height of the cone given by

$$ m^2 = \frac{L \left[ k_1 \cos^2(\psi) + k_3 \sin^2(\psi) - 2k_2 \sin^2(\psi) \cos^2(\psi) - 2k_3 \sin^4(\psi) - k_B/2 \right] - M\sigma}{L \left[ k_2 \sin^2(\psi) \cos^2(\psi) + k_3 \sin^4(\psi) + k_B/4 \right] + M\sigma} \quad (S38) $$
## 2.3: Transitions between deformations

### Transition to conical surface in single elastic constant limit

In the single constant limit, we set $k_1 = k_2 = k_3 = k/3$ and we can re-write Eq. S38 as

$$m^2 = \frac{(4k/3k_B)[1 - 2\sin^2(\psi)] - 2 - 4M\sigma/Lk_B}{(4k/3k_B)\sin^2(\psi) + 1 + 4M\sigma/k_B L} \quad (S39)$$

The results of this are shown in Fig. S1a, which is comparable to Fig 1b in the main text.

If we set $A = 4k/3k_B$ and $B = 4M\sigma/Lk_B$ we arrive at Eq. 3 in the main text. If we consider an aster, $\psi = 0$, and set $m = 0$ we identify the point at which the conical surface becomes stable in the single constant limit, which gives Eq. 4 and the red line on Fig. 1b in the main text.

### Transition to conical surface with variable elastic constants and predetermined phase

For a given set of elastic coefficients and phase we calculate $m^2$ from Eq. S38 (if $m^2 < 0$ we set $m = 0$). This is plotted in Fig. S1b. The green line on Fig. 1c shows where $m^2 = 0$ according to Eq. S38.

### Transition to conical surface with variable elastic constants and variable phase

We assume a conical aster configuration, $\psi = 0$, thus Eq. S38 becomes

$$m^2 = \frac{L[k_1 - k_B/2] - M\sigma}{Lk_B/4 + M\sigma} \quad (S40)$$

The results of this equation are plotted in Fig. S1c. To obtain the transition curves to either an aster or a vortex on a flat disc. We then substitute the solution of Eq. S40 into Eq. S37 to get the energy of a conical surface with an aster.

The conical aster energy is then compared to the energy of an aster or vortex on a flat disc, which are given by:

$$F_{\text{aster,flat}} = Lk_1 + M\sigma \quad (S41)$$

$$F_{\text{vortex,flat}} = Lk_3 + M\sigma \quad (S42)$$

This gives the transition curves in Fig. S1c, which correspond to Eq. 4 and Eq. 5 in the main text, respectively.

<figure><img src="image_6.png" /><figcaption>FIG. S1. Result of Eq. S38 for three sets of variables. (a) Single Frank elastic constant ($k_1 = k_2 = k_3 = k/3$) with and aster ($\psi = 0$) varying $\sigma$, $k$ and $k_B$ under the constraint $k + k_B + \sigma = 1$. This is analogous to Fig. 1b (main text). (b) Varying Frank elastic constants and predetermined phase with $k_B = \sigma = 1/10$, $k_2 = 0$, $k_1 + k_3 = 1$. This is analogous to Fig. 1c (main text). (c) Varying Frank elastic constants with aster defect $\psi = 0$ with $k_B = \sigma = 1/10$, $k_1 + k_2 + k_3 = 1$. This is analogous to Fig. 2a (main text).</figcaption></figure>

## Section 2.4: The limit $\Delta \to 0$

We describe the elastic energy of the nematic field to be given by the Frank free energy, which diverges at the core of the defect. Since we do not model the physics of the core of the defect, this necessitates a finite defect core radius, $\Delta$. Thus, we restrict our view to systems which are much bigger than the defect core radius, i.e. $\lim_{\Delta \to 0} \Delta \to 0$.

In the conical calculation above, Eq. S37, $\Delta$ modulates the relative importance of the surface tension, $\sigma$. We can see this by rescaling the energy by the defect core energy, $L = 2\pi \log(1/\Delta)$. In the limit $\Delta \to 0$, the surface tension becomes insignificant relative to the other contributing energies. This allows us to evaluate all transitions in the $\Delta = 0$ limit by considering the case where $\sigma = 0$. This shows very little change with respect to our previous calculations,
see Fig. S2. There are two noticeable differences: first, the transition between buckled and flat asters in the single Frank constant limit is now independent of the value of $\sigma$. Second, the magnitude of the buckled surfaces is generally higher, this is because there is more energy available to drive the deformation.

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

FIG. S2. Result of Eq. S38 for three sets of variables evaluated for $\Delta = 0$. (a) Single Frank elastic constant ($k_1 = k_2 = k_3 = k/3$) with and aster ($\psi = 0$) varying $\sigma$, $k$ and $k_B$ under the constraint $k + k_B + \sigma = 1$. This is analogous to Fig. 1b (main text). (b) Varying Frank elastic constants and predetermined phase with $k_B = \sigma = 1/10$, $k_2 = 0$, $k_1 + k_3 = 1$. This is analogous to Fig. 1c (main text). (c) Varying Frank elastic constants with aster defect $\psi = 0$ with $k_B = \sigma = 1/10$, $k_1 + k_2 + k_3 = 1$. This is analogous to Fig. 2a (main text).

## Section 2.5: Effect of non-zero $\Delta$

We use a value of $\Delta = 10^{-2}R$ in our numerical calculations, as a balance between accuracy and computational feasibility. To assess the validity of this choice, we recreate the results from Fig. S1 for a large range of $\Delta$.

Fig. S3 shows how the boundaries between buckled and flat sheets varies as $\Delta$ is varied over 7 orders of magnitude. It is clear here that the effect is small as $\Delta$ decreases and that our value of $\Delta = 10^{-2}R$ does not introduce significant in the transitions present.

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

FIG. S3. Position of the boundaries between buckled and flat sheets according to Eq. S37 for values of $\Delta$ that vary over 7 orders of magnitude for three sets of variables. (a) Single Frank elastic constant ($k_1 = k_2 = k_3 = k/3$) with and aster ($\psi = 0$) varying $\sigma$, $k$ and $k_B$ under the constraint $k + k_B + \sigma = 1$. This is analogous to Fig. 1b (main text). (b) Varying Frank elastic constants and predetermined phase with $k_B = \sigma = 1/10$, $k_2 = 0$, $k_1 + k_3 = 1$. This is analogous to Fig. 1c (main text). (c) Varying Frank elastic constants with aster defect $\psi = 0$ with $k_B = \sigma = 1/10$, $k_1 + k_2 + k_3 = 1$. This is analogous to Fig. 2a (main text).
# SECTION 3: MINIMIZING THE ENERGY THROUGH MONTE-CARLO METHOD

## Section 3.1: Evaluating the energy

To solve using Monte-Carlo methods, the equations are recast using a cylindrical Monge gauge. The surface is given by the equations

$$X = r\cos(\theta), \qquad (S43)$$

$$Y = r\sin(\theta), \qquad (S44)$$

$$Z = \zeta(r) \qquad (S45)$$

and the director field is given by the function $\psi(r)$. To fully evaluate all terms of the free-energy density on the surface, we need to know the function $\psi$ and its first derivative and $\zeta$ and its first two derivatives.

We do this by discretizing the functions in the radial direction and set the radius of the patch $R = 1$ in all simulations. The space is divided into $N$ discretized points, and the locations of the center of the first and last bins are at $\epsilon = 1/2N$ and $1 - \epsilon$, respectively. There are two boundaries that we need to accommodate.

### Section 3.1.1: Center of the surface

The center of the surface we assume is given by a spherical cap. At this point the energy associated with the director field diverges, and we consider the spherical cap to include no orientational order. Thus the only terms we need to calculate are those associated with the bending energy and tension of the surface. These are given by

$$F_{tip} = \int k_B \mathcal{H}^2 + \sigma dA \qquad (S46)$$

The surface area of a spherical cap is given by $A = 2\pi L^2(1 - \cos(\alpha))$, where $\alpha$ is the angular size of the spherical cap and $L$ is its radius. The spherical cap extends to a radius $\epsilon$ in the cylindrical radial coordinate, which is the position of the first discretization point, at which point it must smoothly transition into the function $\zeta$. After some trigonometry, we can arrive at

$$F_{tip} = 2\pi k_B \left[ 1 - \frac{1}{\sqrt{\zeta'^2(\epsilon) + 1}} \right] + 2\pi\sigma\epsilon^2 \left[ 1 - \frac{1}{\sqrt{\zeta'^2(\epsilon) + 1}} \right] \frac{(1 + \zeta'^2(\epsilon))}{\zeta'^2(\epsilon)} \qquad (S47)$$

### Section 3.1.2: Edge of the surface

The boundary conditions at the edge of the patch are that $\zeta'(r \ge 1) = M$. This is accommodated within the liquid crystal energy and surface tension by including an imaginary point beyond the boundary with the correct boundary conditions. However, for the curvature terms, which rely specifically on the second derivative of the height function, more care is required. This is to prevent a diverging value of $\zeta''$ below the lengthscale of the discretization of the surface.

We define the second derivative of the height as

$$\zeta''(1) = (M - \zeta'(1 - \epsilon))/\epsilon \qquad (S48)$$

Here $\epsilon = 1/2N$ where $N$ is the number of discretized points in $r$. Thus $\zeta'(1 - \epsilon)$ is the evaluated at the last grid point, which is a distance $\epsilon$ from the boundary.

The area element at this point is given by $J = \sqrt{1 + M^2}$ thus we can write the mean curvature at this point as

$$\mathcal{H} = \frac{M + M^3 + (M - \zeta'(1 - \epsilon))/\epsilon}{2J^3} \qquad (S49)$$

We consider $\mathcal{H}$ to be constant in this small region close to the boundary which has an area equal to $2\pi\epsilon J$. Which results in an additional energy given by

$$F_{edge} = 2\pi k_B \epsilon \frac{(M + M^3 + (M - \zeta'(1 - \epsilon))/\epsilon)^2}{4J^5} \qquad (S50)$$
## Section 3.2: Minimizing the energy through Monte Carlo approach

We assume that the radial height profile of the surface can be approximated by an infinite series of Gaussian functions. This is true for any $L^1$ function (see Yves Myers' "Wavelents and Operators" for a proof). Thus we assume that the function $\zeta$ can be written

$$ \zeta(r) = \sum_{i} A_{i} \exp[-(r - \mu_{i})^{2}/2\sigma_{i}^{2}]. \qquad (S51) $$

Where $A_i$, $\sigma_i^2$ and $\mu_i$ are the amplitude, variance and mean of the constituent Gaussian functions, respectively.

From this we introduce the perturbation

$$ \Delta\zeta(r) = A \exp[-(r - \mu)^{2}/2\sigma^{2}] - A \exp[-(1 - \mu)^{2}/2\sigma^{2}]. \qquad (S52) $$

We have introduced a constant here that fixes $\zeta(1) = 0$.

In the case where the gradient at the boundary is fixed such that $\zeta'(r \ge 1) = M$ we modify the perturbation to satisfy this boundary condition

$$ \Delta\zeta(r) = A \exp[-(r - \mu)^{2}/2\sigma^{2}] + A \exp[-(r - 2 + \mu)^{2}/2\sigma^{2}] - A \exp[-(1 - \mu)^{2}/2\sigma^{2}] - A \exp[-(1 - 2 + \mu)^{2}/2\sigma^{2}]. \qquad (S53) $$

This perturbation has both $\Delta\zeta(1) = 0$ and $\Delta\zeta'(1) = 0$, thus when combined with the initial condition $\zeta(r) = Mr$, the boundary conditions are preserved. It should be noted that this form does not strictly enforce the boundary conditions due to the discretized nature of the computational approach, it is merely used to speed up the Monte Carlo algorithm. The boundary is constrained by the additional edge energy introduced above.

We take a similar approach with the orientation field. We assume the orientation field can be expressed

$$ \psi(r) = \psi_{1} + \sum_{i} B_{i} \exp[-(r - \rho_{i})^{2}/2\Sigma_{i}^{2}]. \qquad (S54) $$

We use the following perturbation $\Delta\psi$

$$ \Delta\psi(r) = B \exp[-(r - \rho)^{2}/2\Sigma^{2}]. \qquad (S55) $$

Similar to previous, when the boundary conditions on the director field are given by $\psi(r \ge 1) = \psi_1$ we use the perturbation

$$ \Delta\psi(r) = B \exp[-(r - \rho)^{2}/2\Sigma^{2}] + B \exp[-(r - 2 + \rho)^{2}/2\Sigma^{2}] - B \exp[-(1 - \rho)^{2}/2\Sigma^{2}] - B \exp[-(1 - 2 + \rho)^{2}/2\Sigma^{2}]. \qquad (S56) $$

This additionally fixes $\Delta\psi(1) = \Delta\psi'(1) = 0$. As before, this is not strictly necessary but is used to increase the speed of the Monte Carlo algorithm.

These perturbations are applied sequentially to the functions $\zeta$ and $\psi$ and are accepted or rejected according to a Boltzmann factor with a gradually reducing temperature in a standard simulated annealing process to find minima in the total free-energy. The parameters are sampled with a uniform distribution with the following limits.

$$ A \in [-0.1, 0.1] \qquad (S57) $$

$$ \mu \in [-1, 1] \qquad (S58) $$

$$ \sigma^2 \in [0.1, 2] \qquad (S59) $$

$$ B \in [-0.1, 0.1] \qquad (S60) $$

$$ \rho \in [-1, 1] \qquad (S61) $$

$$ \Sigma^2 \in [0.1, 2] \qquad (S62) $$

# SECTION 4: DISCUSSION ON BOUNDARY CONDITIONS AND TOPOLOGICAL CONSTRAINTS.

In this work, we studied patches with rotational symmetry, which enforce the existence of a charge +1 topological defect at their core. This is to reflect topological defects that have been associated with morphological changes in
experimental systems [37–40, 54]. Here, we address the question of how these patches can be tiled together to represent a closed surface, such as those often observed in naturally occurring biological systems.

The total topological charge of a vector field on a closed surface is constrained according to the Poincaré-Hopf Theorem, which states that the sum of the topological charges of a vector field is equal to the Euler characteristic of the surface on which they lie. Therefore a simple surface with a spherical topology can be tiled with two +1 defects as the Euler characteristic of a spherical surface is +2. More generally, if we want to tile a spherical surface with many patches containing +1 topological defects, additional negative charges are required to meet the constraints of the Poincaré-Hopf Theorem. This is observed in regenerating freshwater Hydra, on which morphogenesis appears to be driven by +1 defects which are placed according to chemical morphogens, and the remaining -1/2 defects arise as a result of the Poincaré-Hopf Theorem and the nematic nature of the interactions between the filaments [37, 40].

Interestingly, a similar statement can be made about the Gaussian curvature, the integral of which is also proportional to the Euler characteristic of the surface (for a closed surface) via the Gauss-Bonnet Theorem. Thus the total topological charge on a surface and the total Gaussian curvature are equal [55]. Therefore topological defects on closed surfaces merely redistribute the Gaussian curvature that is already present on the surface. In this scenario, the additional negative defects necessitated by the Poincaré-Hopf Theorem can simply be placed in the negative Gaussian curvature regions of the surface necessitated by the Gauss-Bonnet Theorem; which reduces their internal energy. This is again consistent with observations on freshwater Hydra, on which the -1/2 defects occur in the regions of maximally negative Gaussian curvature [29, 37, 40].

## SECTION 5: PARAMETERS AND BOUNDARY CONDITIONS FOR ALL NUMERICAL RESULTS

All simulations are performed on a discretized line of $N = 100$ points. Parameter values and boundary conditions for all numerical work is provided in the tables below.

<table><thead><tr><th>Figure</th><th>Parameters</th><th>Constraints</th></tr></thead><tbody><tr><td>Fig. 1b</td><td>$k_1 = k_2 = k_3 = k/3$<br>$\sigma + k + k_B = 1$</td><td>$\zeta''(R) = 0$<br>$\psi''(R) = 0$</td></tr><tr><td>Fig. 1c</td><td>$k_1 + k_3 = 1$<br>$\sigma = k_B = 1/10$<br>$k_2 = 0$</td><td>$\zeta''(R) = 0$<br>$\psi'(r) = 0$</td></tr><tr><td>Fig. 1d</td><td>$k_1 = k_3 = 1/2$<br>$\sigma = k_B = 1/10$<br>$k_2 = 0$</td><td>$\zeta''(R) = 0$<br>$\psi(r) = 0$<br>$\psi'(r) = 0$</td></tr><tr><td>Fig. 1e</td><td>$k_3 = 1$<br>$\sigma = k_B = 1/10$<br>$k_1 = k_2 = 0$</td><td>$\zeta''(R) = 0$<br>$\psi(r) = \pi/8$<br>$\psi'(r) = 0$</td></tr></tbody></table>

<table><thead><tr><th>Figure</th><th>Parameters</th><th>Constraints</th></tr></thead><tbody><tr><td>Fig. 2a</td><td>$k_1 + k_2 + k_3 = 1$</td><td>$\zeta''(R) = 0$</td></tr><tr><td>Fig. 2b</td><td>$\sigma = k_B = 1/10$</td><td>$\psi''(R) = 0$</td></tr><tr><td>Fig. 2c</td><td>$k_1 = 1/2$<br>$k_2 = 1/10$<br>$k_3 = 4/10$<br>$\sigma = k_B = 1/10$</td><td>$\zeta''(R) = 0$<br>$\psi''(R) = 0$</td></tr><tr><td>Fig. 2d</td><td>$k_1 = 13/20$<br>$k_2 = 1/10$<br>$k_3 = 1/4$<br>$\sigma = k_B = 1/10$</td><td>$\zeta''(R) = 0$<br>$\psi''(R) = 0$</td></tr></tbody></table>
<table><thead><tr><th>Figure</th><th>Parameters</th><th>Constraints</th></tr></thead><tbody><tr><td>Fig. 3a</td><td>$k_1 + k_2 + k_3 = 1$</td><td>$\zeta'(R) = 0$</td></tr><tr><td>Fig. 3b</td><td>$k_2 = 1/3$<br>$\sigma = k_B = 1/10$</td><td>$\psi(R) = \psi_R$<br>$\psi'(R) = 0$</td></tr></tbody></table>

<table><tbody><tr><td>Fig. 3c</td><td>$k_1 = 1/2$</td><td>$\zeta'(R) = 0$</td></tr><tr><td>Fig. 4b</td><td>$k_2 = 1/3$<br>$k_3 = 1/6$<br>$\sigma = k_B = 1/10$</td><td>$\psi(R) = 0$<br>$\psi'(R) = 0$</td></tr></tbody></table>

<table><tbody><tr><td>Fig. 3d</td><td>$k_1 = 2/3$</td><td>$\zeta'(R) = 0$</td></tr><tr><td>Fig. 4c</td><td>$k_2 = 1/3$<br>$k_3 = 0$<br>$\sigma = k_B = 1/10$</td><td>$\psi(R) = 0$<br>$\psi'(R) = 0$</td></tr></tbody></table>

<table><tbody><tr><td>Fig. 3e</td><td>$k_1 = 0$<br>$k_2 = 1/3$<br>$k_3 = 2/3$<br>$\sigma = k_B = 1/10$</td><td>$\zeta'(R) = 0$<br>$\psi(R) = \pi/2$<br>$\psi'(R) = 0$</td></tr></tbody></table>

<table><thead><tr><th>Figure</th><th>Parameters</th><th>Constraints</th></tr></thead><tbody><tr><td>Fig. 4a</td><td>$k_1 = 1/2$</td><td>$\zeta'(R) = 0$</td></tr><tr><td>Black</td><td>$k_3 = 1/6$<br>$\sigma = k_B = 1/10$</td><td>$\psi(R) = 0$<br>$\psi'(R) = 0$</td></tr></tbody></table>

<table><tbody><tr><td>Fig. 4a</td><td>$k_1 = 2/3$</td><td>$\zeta'(R) = 0$</td></tr><tr><td>Cyan</td><td>$k_3 = 0$<br>$\sigma = k_B = 1/10$</td><td>$\psi(R) = 0$<br>$\psi'(R) = 0$</td></tr></tbody></table>

<table><tbody><tr><td>Fig. 4a</td><td>$k_1 = 0$</td><td>$\zeta'(R) = 0$</td></tr><tr><td>Magenta</td><td>$k_3 = 2/3$<br>$\sigma = k_B = 1/10$</td><td>$\psi(R) = \pi/2$<br>$\psi'(R) = 0$</td></tr></tbody></table>

<table><thead><tr><th>Figure</th><th>Parameters</th><th>Constraints</th></tr></thead><tbody><tr><td>Fig. 4d</td><td>$k_1 = 1/2$<br>$k_2 = 1/3$<br>$k_3 = 1/6$<br>$\sigma = k_B = 1/10$</td><td>$\zeta'(R) = 1$<br>$\psi(R) = 0$<br>$\psi'(R) = 0$</td></tr></tbody></table>

<table><tbody><tr><td>Fig. 4e</td><td>$k_1 = 2/3$<br>$k_2 = 1/3$<br>$k_3 = 0$<br>$\sigma = k_B = 1/10$</td><td>$\zeta'(R) = 1$<br>$\psi(R) = 0$<br>$\psi'(R) = 0$</td></tr></tbody></table>

<table><tbody><tr><td>Fig. 4f</td><td>$k_1 = 0$<br>$k_2 = 1/3$<br>$k_3 = 2/3$<br>$\sigma = k_B = 1/10$</td><td>$\zeta'(R) = 1$<br>$\psi(R) = \pi/2$<br>$\psi'(R) = 0$</td></tr></tbody></table>

<table><thead><tr><th>Figure</th><th>Parameters</th></tr></thead><tbody><tr><td>Fig. S1a</td><td>$k_1 = k_2 = k_3 = k/3$<br>$k + \sigma + k_B = 1$<br>$\psi = 0$<br>$\Delta = 1/100$</td></tr></tbody></table>

<table><tbody><tr><td>Fig. S1b</td><td>$k_1 + k_3 = 1$<br>$k_2 = 0$<br>$\sigma = k_B = 1/10$<br>$\Delta = 1/100$</td></tr></tbody></table>

<table><tbody><tr><td>Fig. S1c</td><td>$k_1 + k_2 + k_3 = 1$<br>$\sigma = k_B = 1/10$<br>$\psi = 0$<br>$\Delta = 1/100$</td></tr></tbody></table>
SECTION 6: QUANTITATIVE COMPARISON TO BIOLOGICAL SYSTEMS

In this section, we summarise the main arguments to make a comparison with biological systems. The values of the material and geometrical parameters can be found in Table S1. Due to the complexity of the systems considered and the simplifications required for our theoretical treatment, we consider only orders of magnitude.

Several biological systems share the main characteristics of a fluid membrane with nematic order. In particular, a quasi-two dimensional geometry, a fluid-like behaviour in the experimental time scale, and nematic orientational order over macroscopic time and length scales. In the following, we focus on two special cases: films of cytoskeletal filaments and thin tissues.

First, let us consider a thin layer made of cytoskeletal filaments, like actin or microtubules, on top of a supported lipid bilayer, [47–49]. For a lipid bilayer, the bending rigidity ranges from $k_B = 10^{-8} – 10^{-7} \mu N\mu m$, and the surface tension ranges from $\sigma = 10^{-6} – 10^{-3} N/m$ [56–59]. Actin and microtubule films exhibit nematic phases [47–49], and the Frank constant $K$ for both filament types has been measured. Specifically, a thin film of actin filaments on an oil-water interface has a Frank constant ranging from $K = 10^{-1} – 1 pN$, [50, 51]. The reduced Frank constant $k = Kh$ depends also on the thickness of the layer $h$. Taking $h = 0.1 \mu m$, which is typical thickness of an actin cortex in cells [60–62], leads to $k = 10^{-8} – 10^{-7} \mu N\mu m$. For thin films of microtubules on an oil-water interface, the reduced Frank constant was $k = 10^{-3} \mu N\mu m$ and the layer thickness was $h = 10 \mu m$ [52]. The corresponding Frank constant for the latter case is compatible with that found in spindles of Xenopus egg [63].

Next, let us consider a monolayer (or thin multilayer) of cells that can be supported by a thin elastic substrate, [37, 64, 65]. For epithelial cells, the bending rigidity is measured as $k_B = 10^{-1} – 1 \mu N\mu m$, and the tension can vary from $\sigma = 10^{-3} – 10^{-1} N/m$ [64, 66–69]. Several types of cells, including epithelial or fibroblast, can exhibit nematic phases, [70–72]. At time of writing, neither the Frank constant nor the reduced Frank constant has been directly measured in tissues, but several works estimated their values from indirect experimental measurements. The value of the Frank constant depends on the cell type, being on the range of $K = 10^3 – 10^4 pN$ for epithelial cells [73, 74] and on the scale of $K = 10^5 pN$ for fibroblast cells [46, 75]. Taking $h = 10 \mu m$, which is the typical thickness of a cell monolayer, leads to the values of the reduce Frank constant on the range of $k = 10^{-2} – 10^{-1} \mu N\mu m$ for epithelial cells and on the scale of $K = 1 \mu N\mu m$ for fibroblast cells.

We show that the instability threshold of a flat disc with an aster at its center is found to depend on the elastic parameters of the fluid membrane and embedded nematic field, in particular the surface tension $\sigma$, the bending rigidity $k_B$, and the reduced Frank constant $k$, as well as, geometrical parameters, such as the radius of the disc $R$. For simplicity, we consider that all the reduced Frank constants are of the same order of magnitude, thus we take the one-constant approximation (i.e. $k_1 = k_2 = k_3 \sim k$). Next, we evaluate this threshold for the previous biological systems in two limiting cases controlled by the ratio $k_B/\sigma R^2$.

In the bending dominated regime, $k_B \gg \sigma R^2$, the instability threshold is dependent on the ratio $k/k_B$ and the critical magnitude is of order of $(k/k_B)_c \sim 1$, see Eq. (4) in the main text. For a film of actin filaments, we expect that this ratio is $k/k_B = 1$. If the film is made of microtubules, then the ratio ranges from $k/k_B = 10^4 – 10^5$, which

<table><caption>TABLE S1. Table of material parameters and geometrical parameters of some biological systems. The symbol * denotes the values that were estimated. For a thin layer, the reduced Frank constant $k = K * h$ depends on the Frank constant $K$ and the thickness of the layer $h$.</caption><thead><tr><th></th><th>Surface tension<br>σ[N/m]</th><th>Bending rigidity<br>k<sub>B</sub>[μNμm]</th><th>Length Scale<br>√k<sub>B</sub>/σ[μm]</th><th>Reference</th></tr></thead><tbody><tr><td>Lipid vesicles</td><td>10<sup>-6</sup> – 10<sup>-3</sup></td><td>10<sup>-8</sup> – 10<sup>-7</sup></td><td>10<sup>-2</sup> – 10<sup>-1</sup> (*)</td><td>[56–59]</td></tr><tr><td>Epithelial cell monolayer</td><td>10<sup>-3</sup> – 10<sup>-1</sup></td><td>10<sup>-1</sup> – 1</td><td>1 – 10 (*)</td><td>[64, 66–69]</td></tr><tr><td></td><td>Frank constant<br>K[pN]</td><td>Thickness<br>h[μm]</td><td>Reduced Frank constant<br>k[μNμm]</td><td>Reference</td></tr><tr><td>F-actin film<br>at the oil-water interface</td><td>10<sup>-1</sup> – 1</td><td>10<sup>-1</sup></td><td>10<sup>-8</sup> – 10<sup>-7</sup> (*)</td><td>[50, 51]</td></tr><tr><td>Microtubule film<br>at the oil-water interface</td><td>10<sup>2</sup> (*)</td><td>10</td><td>10<sup>-3</sup></td><td>[52, 76]</td></tr><tr><td>Xenopus egg<br>extract spindles</td><td>10<sup>2</sup></td><td>-</td><td>-</td><td>[63]</td></tr><tr><td>Epithelial cell monolayer</td><td>10<sup>3</sup> – 10<sup>4</sup> (*)</td><td>10</td><td>10<sup>-2</sup> – 10<sup>-1</sup> (*)</td><td>[73, 74]</td></tr><tr><td>Fibroblast cell monolayer</td><td>10<sup>5</sup> (*)</td><td>10</td><td>1 (*)</td><td>[46, 75]</td></tr></tbody></table>
is above the instability threshold. In microtubule films the thickness of the layer can be modulated and may allow to control the value of this ratio $k/k_B$ [52, 76]. For tissues and ignoring the elasticity of the underlying substrate, the ratio is expected to be $k/k_B = 10^{-1} - 1$.

In the tension dominated regime, $k_B \ll \sigma R^2$, the instability threshold of a flat interface is controlled by the ratio $k/\sigma R^2$, and the critical value is of order $(k/\sigma R^2)_c \sim 1$. One can estimate an upper bound by replacing the radius $R$ with the layer thickness $h$, which is the smallest lengthscale in the problem. Depending on the surface tension $\sigma$, the upper bound is $k/\sigma h^2 = 10^{-3} - 10$ for a film of actin filaments and $k/\sigma h^2 = 10^{-2} - 10$ for a film of microtubules. Ignoring the elasticity of the underlying substrate, the upper bound for a cell monolayer is $k/\sigma h^2 = 10^{-3} - 10$. In all cases, the upper bound is at most a factor 10 larger than the critical value, which suggest that in low-tension regimes, these cases can be close to the instability threshold.

Note that because the geometry is quasi-two dimensional, the thickness of the layer $h$ needs to be smaller than the typical radius $R$. Therefore we can estimate an upper bound for the ratio $k_B/\sigma R^2$. For both biological examples, we found that $k_B/\sigma h^2 \le 1$, suggesting that experimental cases are in a tension dominated regime $k_B/\sigma R^2 \ll 1$.

To summarise, in this section we evaluated the proximity to the instability threshold of a flat fluid membrane with an aster topological defect for two biological systems: films of cytoskeletal filaments on a supported lipid bilayer and cell monolayers on a supported elastic substrate. In the bending dominated regime $k_B \gg \sigma R^2$, both films of actin filaments and cell monolayers can be close to the instability threshold, and films of microtubules are above the instability threshold. In the tension dominated regime $k_B \ll \sigma R^2$, we found that all these biological examples can be close to the instability threshold. Our arguments here focus on orders of magnitude, to know whether a flat membrane can be rendered unstable for the cases that are close to the instability threshold would depend on the precise values of the elastic and geometrical parameters. Two possible ways to cross the instability threshold in these cases is to either operate in the bending dominated regime, by for instance reducing the surface tension of the supporting substrates, or to increase the reduced Frank constant by for instance increasing the number density of cytoskeletal filaments or cells.