Curvature controlled defect dynamics in topological active nematics

TL;DR

This paper investigates how topological constraints and geometry influence defect dynamics in active nematic films, revealing tunable oscillatory and rotating states linked to curvature and shape.

Contribution

It introduces a detailed analysis of defect behavior in topologically constrained active nematics on ellipsoids, highlighting the role of geometry in defect dynamics.

Findings

Defects exhibit two main dynamic modes: oscillatory and rotating.

Defect locations correlate with high Gaussian curvature and umbilical points.

Limits are identified for modeling defects as self-propelled particles.

Abstract

We study the spatiotemporal patterns that emerge when an active nematic film is topologically constraint. These topological constraints allow to control the non-equilibrium dynamics of the active system. We consider ellipsoidal shapes for which the resulting defects are 1/2 disclinations and analyze the relation between their location and dynamics and local geometric properties of the ellipsoid. We highlight two dynamic modes: a tunable periodic state that oscillates between two defect configurations on a spherical shape and a tunable rotating state for oblate spheroids. We further demonstrate the relation between defects and high Gaussian curvature and umbilical points and point out limits for a coarse-grained description of defects as self-propelled particles.

Figures

Equations

$$\eta =\frac{1}{4Nc}\sum _{i=1}^4\sum _{t=ts}^{te}|h_i(t)|,$$

$$\begin{array}[]{ll}\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{q}=\mathbf{v},\quad\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{v}=\mathbf{F}-G\left(\mathbf{q}\right)^{T}\bm{\lambda},\quad\displaystyle g\left(\mathbf{q}\right)=\mathbf{0}\end{array}$$

$${\mathbf{F}}_i=-\gamma {\mathbf{v}}_i+\sum _{j=1}^N{\mathbf{F}}_{ij}+{\mathbf{F}}_i^{ac}$$

$$\begin{array}[]{ll}\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{n}_{i}=\bm{\omega}_{i}\times\mathbf{n}_{i},\quad\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\bm{\omega}_{i}=-\gamma_{a}\bm{\omega}_{i}+\mathbf{T}_{i}\left(\mathbf{q},\mathbf{n}\right)\end{array}$$

$${\mathbf{F}}_i^{ac}=-v_0\sum _{j\in U(i)}{\mathbf{Q}}^j\frac{{\mathbf{q}}_i-{\mathbf{q}}_j}{q_{ij}^2}.$$

$${\mathbf{T}}_i=J\sum _{j\in U(i)}\left(\left({\mathbf{n}}_i\cdot {\mathbf{n}}_j\right)\left({\mathbf{n}}_i\times {\mathbf{n}}_j\right)\right).$$

$${\mathbf{T}}_i=J\sum _{j\in U(i)}\left(1+{\mathbf{n}}_i\cdot {\mathbf{n}}_j\right)\frac{{\mathbf{n}}_i\times {\mathbf{n}}_j}{|{\mathbf{n}}_i\times {\mathbf{n}}_j{|}^2}$$

$$K=\frac{a^2{b}^6{c}^6}{{\left(c^4{b}^4+c^4(a^2-b^2)y^2+b^4(a^2-c^2)z^2\right)}^2}.$$

$${\left(\pm a\sqrt{\frac{a^2-b^2}{a^2-c^2}},0,\pm c\sqrt{\frac{b^2-c^2}{a^2-c^2}}\right)}^T$$

$$P_i=\frac{1}{{\sum }_j{w}_{ij}}\sum _j\frac{w_{ij}}{2}(3{\mathbf{n}}_i\cdot {\mathbf{n}}_j-1)$$

Full text

Active systems are characterized by constant input of energy, which is converted by autonomous constituents into directed motion, leading to spatiotemporal patterns. These phenomena range from the macroscale, e.g. flocks of birds Cavagnaetal_PNAS_2010 or schools of fish Hemelrijketal_BE_2005 to the microscale, e.g. bacterial colonies Wensinketal_PNAS_2012 , migrating tissue cells Szaboetal_PRE_2006 or active nematic films Giomietal_PRL_2013 . If such systems are confined on curved surfaces, topological constraints strongly influence the emerging spatiotemporal patterns. Using these topological constraints to guide collective cell behavior might be a key in morphogenesis and active nematic films on surfaces have been proposed as a promising road to engineer synthetic materials that mimic living organisms Keberetal_Science_2014 . However, the complex dynamics of such topological active systems remains wildly unexplored. As in passive systems the mathematical Poincaré-Hopf theorem forces topological defects to be present in the nematic film. On a sphere this leads to an equilibrium defect configuration with four +1/2 disclinations arranged as a tetrahedron Lubenskietal_JP_1992 ; Shinetal_PRL_2008 ; Skacejetal_PRL_2008 , see Figure 1 The disclinations repel each other and this arrangement maximizes their distance. In active systems unbalanced stresses drive this configuration out of equilibrium. But in contrast to planar active nematics with continuous creation and annihilation of defects Sanchezetal_Nature_2012 ; DeChampetal_NatureM_2015  the creation of additional defect pairs can be suppressed on curved surfaces, which is demonstrated in Keberetal_Science_2014 for an active nematic film of microtubules and molecular motors, encapsulated within a spherical lipid vesicle. This provides an unique way to study the dynamics of the four defects in a controlled manner and led to the discovery of a tunable periodic state that oscillates between the tetrahedral and a planar defect configuration. We confirm this finding by computer simulations, see Figure 1.

Within a coarse-grained model +1/2 disclinations in planar active nematic films can be effectively described by self-propelled particles with a velocity proportional to the activity Giomietal_PRL_2013 . In Keberetal_Science_2014 this relation is extended to spherical nematics. Four self-propelled particles on a sphere also oscillate between the planar and tetrahedral configuration. Both descriptions can be quantitatively linked to each other, but also differences can be pointed out, which become more evident for more general surfaces. For non-constant Gaussian curvature constraints local geometric properties influence the position of the defects and thus can be used to control defect dynamics. We are concerned with a systematic investigation of the impact of such constraints on the emergence of complex patterns and oscillations.

Results

For active systems in flat geometries various theoretical descriptions have been proposed, see e.g. Marchettietal_RMP_2013 ; Bechingeretal_RMP_2016 . One of the most studied approaches are Vicsek-like models Visceketal_PRL_1995 . They consider particles, which travel at a constant speed to represent self-propulsion, whose direction changes according to interaction rules which comprise explicit alignment and noise. In contrast to equilibrium systems long-range order emerges for two dimensional systems with low noise. We consider an extension of these models which includes excluded volume Tailleuretal_PRL_2008 ; Henkesetal_PRE_2011 ; Filyetal_PRL_2012 ; Bialkeetal_PRL_2012 and classify systems by the head-tail symmetry of their particles in polar or nematic. For active polar particles these models have been formulated on a sphere Sknepneketal_PRE_2015 and on ellipsoidal surfaces Ehrigetal_arXiv_2016 . In these situations a robust rotating-band structure around the waist, with two +1 defects at the poles is found on a sphere. On an ellipsoid the location of the defects is linked to local geometric properties, similar to vortices in surface fluids Turneretal_RMP_2010 ; Nitschkeetal_JFM_2012 ; Reutheretal_MMS_2015 ; Nitschkeetal_arXiv_2016 . The defects are related to the Gaussian curvature and to the umbilical points of the surface (see Materials and Methods for a geometric description). For spheroidal ellipsoids there are two umbilical points, which locate the two +1 defects. This configuration is more stable for prolate spheroids, where the umbilical points are at the points of maximal Gaussian curvature at the poles and less stable for oblate spheroids, where the umbilical points and the maximum in Gaussian curvature are separated. As in the spherical case a rotating-band structure is formed, with possible sub-bands which counter rotate depending on the initial condition. New dynamical features are found for non-spherical ellipsoids. They have four umbilical points. For lower velocities the defects encircle pairs of umbilical points and for larger velocities the defects are found at the high Gaussian curvature regions between each pair of umbilical points. With this richness in dynamics found for active polar particles on non-constant Gaussian curvature surfaces, we expect similar behavior for active nematic particles and ask up to which complexity of the geometry the dynamics of the four 1/2 disclinations can be effectively described by self-propelled polar particles.

To answer these questions, we first analyze the spherical case in more detail. In addition to the oscillation between the planar and tetrahedral defect configuration on a spherical vesicles and a tunable frequency by the activity and self-propulsion velocity we also track the positions of the defects. Computing the power spectrum from the time series for the average angle $<\alpha>=\frac{1}{6}\sum_{i<j}\alpha_{ij}$ we obtain the frequency for the oscillations, which linearly depend on the activity. The same results, but with a small offset and a different slope are obtained for the coarse-grained description by self-propelled particles, see Figure 2. As a consequence for each activity in the nematic film a self-propulsion velocity can be determined in the coarse-grained description, which resamples the frequency of the planar-tetrahedral defect oscillation. Differences between both descriptions are found if we compare the trajectories of the defects and self-propelled particles. Within the considered time interval the 1/2 disclinations are locally confined, each defect only covers part of the vesicle. This is in contrast to the trajectories of the self-propelled particles, which rotate within a band structure leaving parts of the vesicle uncovered, see Figure 2. The experimental defect trajectories in Keberetal_Science_2014 differ from both descriptions, they are global, covering the whole vesicle. The discrepancy might be a consequence of the considered short-range interactions in the model for the active nematic film.

We next consider spheroidal ellipsoids. They are characterized by the aspect ratio $a/c$ and $a=b$, with $a$, $b$ and $c$ the length of the major axis. Due to the symmetry all geometric properties can be characterized with respect to the polar axis. As the geometry is topologically equivalent to a sphere we expect for passive systems again a minimal energy configuration with four 1/2 disclinations. They still try to maximize their distance, but are now also influenced by local geometric properties. The 1/2 disclinations tend to accumulate in regions of high Gaussian curvature Batesetal_SM_2009 ; Serra_LC_2016 . Computer simulations for thin nematic shells have shown that for prolate ellipsoids pairs of defects are located at opposite ends close to the poles. The defects in each pair arrange at opposite sides of the surface and tend to align perpendicular to the pair at the other pole Batesetal_SM_2009 . As the distance between the defects is no longer maximized, the geometric effect seems to dominate the repulsion in this case. For oblate ellipsoids the 1/2 disclinations are found near the waist, where the Gaussian curvature is largest. Again two pairs of defects are found, one on each side. They repel each other and are mutually perpendicular to the other pair, leading to an alternating ring of 1/2 disclinations, one above and one below the waist. This behavior seems to be independent of the film thickness Batesetal_SM_2009 , we have confirmed this behavior by our surface model without activity.

For active systems we observe again oscillatory behavior, see Figure 3. For prolate spheroids ($a/c<1$) only two 1/2 disclinations are located at the poles, whereas the other two oscillate around the waist. The oscillations are very noisy and can not be tuned by the activity. Even if the distance between the two 1/2 disclinations at the waist is not optimal the average distance between all four defects is larger than in the passive case. While the 1/2 disclinations are still attracted by the high curvature regions at the poles, the active forces push one of the defects away leading to the observed metastable configuration. Within a transition zone ($a/c\approx 1$) we observe similar behavior as in the spherical case ($a/c=1$) without any defect localization. The behavior changes for oblate spheroids ($a/c>1$), where all four 1/2 disclinations are along the waist, maintaining a maximal distance to each other. This behavior is similar to the passive system. However, the defects now oscillate between both sides. The frequency of the alternating oscillations above and below the waist can be extracted for various activities. However, a clear functional dependency on the activity could not be found. If the aspect ratio is further increased the situation changes to pairs of 1/2 disclinations which rotate around the umbilical points at the poles. The defects are no longer located at positions of maximal Gaussian curvature. The high curvature value at the waist creates a distortion of the nematic film, which can be seen from the nematic order parameter. It somehow serves as a barrier for the 1/2 disclinations preventing them from crossing the waist. The rotation is a consequence of the activity and the unfavorable short distance with respect to each other. The frequency of the rotation depends on the activity and can be tuned, see Figure 4. Also the transition to this rotating state depends on the strength of the activity. As stronger the activity as longer it is possible for the defects to cross the barrier at the waist. A tendency to locate the defects away from the high Gaussian curvature waist can also be seen for the passive case.

The four different regimes are shown in Figure 5. using the order parameter

[TABLE]

with $N$ the number of particles, $[ts,te]$ an appropriate time interval and $h_{i}(t)$ the height of the defect $i$ along the polar axis with respect to the waist at time $t$. We have $\eta=1$ if all defects are at the poles, $\eta=0$ if they are at the waist and $\eta=0.5$ if they are homogeneously distributed along the polar axis.

Within the coarse-grained description by self-propelled polar particles, using the corresponding self-propulsion velocity according to Figure 2, we obtain a qualitatively different behavior. Within the considered parameter regime, the values for $\eta$ are independent of the self-propulsion velocity. For aspect ratios $a/c<0.5$ the particles rotate on closed trajectories, well separated from each other at approximately equal distance along the polar axis. The transition zone with sphere-like behavior is more extended than for the nematic defects. For $0.5<a/c<2$ a band structure is formed around the waist, which shrinks with increasing aspect ratio. For $a/c>2$ all particles are positioned at the waist, rotating in one direction and maintaining their distance. The regime with pairwise rotating defects around the umbilical points could not be found within the coarse-grained model (A more detailed description is given in the SI).

Non-spherical ellipsoids, which are characterized by $a\neq b$, $a\neq c$ and $b\neq c$, have four umbilical points. They are either prolate-like or oblate-like but in any case have two distinct points of maximal Gaussian curvature. We thus analyze the distance of the four 1/2 disclinations with respect to the umbilical points and the points of maximal Gaussian curvature using the average geodesic distances $<DDU>$ and $<DDG>$, respectively. Figure 6, which is inspired by Ehrigetal_arXiv_2016 shows the distances as a function of the aspect ratios $a/b$ and $a/c$. Spheroids are also included, the first column shows the previous results for oblate and the diagonal for prolate geometries. Each row in between thus corresponds to a transition from oblate-like to prolate-like geometries. In most cases the 1/2 disclinations are closer to the high Gaussian curvature points than to the umbilical points, with the only exception for oblate-like ellipsoids with a large aspect ration $a/c\geq 4$. This leads to the conclusion that 1/2 disclinations tend to be attracted by points of high Gaussian curvature.

Discussions

In Keberetal_Science_2014 it was shown that in a confined active system, a dense suspension of microtubules and molecular motors on the surface of a spherical lipid vesicle, cyclic oscillations between defect configurations can be observed. They result from topological constraints and the coupling between velocity fields and defect-defect interactions. These findings may push forward the design of systems that harness the ability of nanoscale active matter to transform chemical energy into mechanical work. On non-spherical surfaces defects are known to be strongly influenced by local geometric properties. The induced geometric interaction can lead to locating of defects, which is established for vortices in surface fluids Turneretal_RMP_2010 ; Nitschkeetal_JFM_2012 ; Nitschkeetal_arXiv_2016 and vortices, sources and sinks in polar systems Nestleretal_arXiv_2016 ; Ehrigetal_arXiv_2016 . For strong variations in geometric properties it has even be found computationally that lower energy minima in passive systems can be formed by creating additional defects Reutheretal_MMS_2015 ; Nestleretal_arXiv_2016 for surface fluids and surface polar particles, respectively. Our work extends the understanding of the delicate relations between topology, geometry and defect dynamics on non-spherical shapes for the system considered in Keberetal_Science_2014 . We are concerned with ellipsoidal surfaces and identify crucial geometric features which influence collective motion patterns in active nematic films. We have shown that 1/2 disclinations are related to both, maxima in the Gaussian curvature and umbilical points of the surface. On prolate spheroids maxima in Gaussian curvature and umbilical points coincide, they are located at the two poles and attract the 1/2 disclinations. However, the repulsive defect-defect interaction allows only two of the defects to be located at the poles, the other two try to maximize their distance and are located around the waist, where they oscillate. Spherical like shapes lead to similar behavior as observed on a sphere, with no distinguished location of the defects and an oscillation between a tetrahedral and planar defect configuration. For oblate spheroids all 1/2 disclinations are located at the waist, the region of high Gaussian curvature. They again maximize their distance and oscillate. With increasing aspect ratio $a/c$ the situation changes. The defects can no longer cross the waist, where the high Gaussian curvature leads to a distortion of the nematic order. As a consequence pairs of 1/2 disclinations rotate around the umbilical points. The frequency of the rotation depends on the activity and can be tuned. This found rotating state is an other step towards a controllable transformation of chemical into mechanical energy in nanoscale active matter and asks for experimental validation. The results for non-spheroidal ellipsoids confirm these findings, even if the separation of the different states is not as distinct as in Figure 5. A smooth transition of the dynamics between prolate-like and oblate-like shapes is identified in Figure, with a clear tendency of the 1/2 disclinations to locate at points of maximal Gaussian curvature. Only for extrem values of $a/c$ and almost spheroidal shapes the situation changes and the rotating state around the umbilical points could be identified.

We further demonstrate that the proposed coarse-grained description of 1/2 disclinations in active nematic matter by self-propelled particles fails if geometric properties come into play. Already on spherical shapes the trajectories of the defects and the self-propelled particles differ significantly and on spheroidal ellipsoids both descriptions not even qualitatively agree.

In summary we explored the complex interaction of topology, geometry and defect dynamics in nematic films on ellipsoidal surfaces and demonstrated how topological constraints and geometric properties can be used to control the collective behavior in nanoscale active matter. The non-linear coupling between non-constant Gaussian curvature and defect-defect interactions leads to tunable spatiotemporal patterns. Among these findings is a stable rotating state on strongly oblate-like ellipsoids, which suggests an other pathway towards a controllable generation of mechanical work in nanoscale active matter. The richness of physics observed in our work will further increase if the underlying shape is deformable. First experimental results of such an interplay between activity-driven defect motion and deformability of the vesicle are already shown in Keberetal_Science_2014 and discussed in Yeomans_NM_2015 . However, for theoretical descriptions of these phenomena new methods will be required.

Materials and Methods

We consider a more general approach than the Viscek-like models confined on a sphere or an ellipsoids in Sknepneketal_PRE_2015 ; Ehrigetal_arXiv_2016 .

.1 Equations of motion

We consider $N$ active particles of mass $m_{i}=1$, which are constrained to move on a surface algebraically described by $g(\mathbf{q})=0$, with particle positions $\mathbf{q}=\left(\mathbf{q}_{1},\ldots,\mathbf{q}_{N}\right)$. Newton’s equations of motion (EOM) with holomonic constraint $g(\mathbf{q})$ read:

[TABLE]

with forces $\mathbf{F}=(\mathbf{F}_{1},\ldots,\mathbf{F}_{N})$ and velocities $\mathbf{v}=(\mathbf{v}_{1},\ldots,\mathbf{v}_{N})$. $\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{N})$ are the Lagrange multipliers and $G(\mathbf{q})=\nabla_{\mathbf{q}}g(\mathbf{q})$ is the Jacobian of $g(\mathbf{q})$. The force $\mathbf{F}_{i}$ can be written as:

[TABLE]

where $\gamma$ is the translational friction coefficient, $\mathbf{F}_{i}^{ac}$ the active force acting on particle $i$ and $\mathbf{F}_{ij}$ the pair-interaction force between particle $i$ and particle $j$. Additionally every particle has an internal degree of freedom, its orientation $\mathbf{n}_{i}$. Denoting by $\bm{\omega}_{i}$ the angular velocity we have the following EOM for the orientational dynamics:

[TABLE]

where $\gamma_{a}$ is the rotational friction coefficient and $\mathbf{T}_{i}(\mathbf{q},\mathbf{n})$ is the torque acting on particle $i$, with $\mathbf{n} =(\mathbf{n}_{1},\ldots,\mathbf{n}_{N})$. Depending on the specific form for the active force $\mathbf{F}_{i}^{ac}$, the pair-interaction force $\mathbf{F}_{ij}$, the torque $\mathbf{T}_{i}(\mathbf{q},\mathbf{n})$ and the holomonic constraint $g(\mathbf{q})$ we will be able to describe polar and nematic active systems on various surfaces.

.2 Active polar particles

For active polar particles on a sphere of radius $R$ we would specify Sknepneketal_PRE_2015 $\mathbf{F}_{i}^{ac}=v_{0}\mathbf{n}_{i}$ with a constant self-propulsion velocity $v_{0}$, $\mathbf{F}_{ij}=k(2\sigma-q^{g}_{ij})\frac{\mathbf{q}_{i}-\mathbf{q}_{j}}{q_{ij}}$ for $q^{g}_{ij} <2\sigma$ and $\mathbf{F}_{ij}=0$ otherwise, a short-range repulsion between spheres of radius $\sigma$, with elastic constant $k$, euclidian distance $q_{ij}=|\mathbf{q}_{i}-\mathbf{q}_{j}|$ and geodesic distance $q^{g}_{ij}=|\mathbf{q}_{i}-\mathbf{q}_{j}|_{g}$. Parallel orientations between neighbouring particles are favored and therefore we use the aligning torque $\mathbf{T}_{i}(\mathbf{q},\mathbf{n})=-J\sum_{j\in U(i)}\left(\mathbf{n}_{i}\times\mathbf{n}_{j}\right)$, with $J>0$ the strength and $U(i)$ the first shell of neighbors of particle $i$, identified as all the particles within a cutoff radius of $2.4\sigma$ from $\mathbf{r}_{i}$. The holomonic constraint for a sphere of radius $R$ reads $g(\mathbf{q}_{i})=q_{i,1}^{2}+q_{i,2}^{2}+q_{i,3}^{2}-R^{2}$, with $\mathbf{q}_{i}=(q_{i,1},q_{i,2},q_{i,3})\in{\cal R}^{3}$. This approach can be used to reproduce the results in Sknepneketal_PRE_2015 in which the overdamped limit, the euclidian distance instead of the geodesic distance and an additional noise term are considered.

.3 Active nematic particles

To describe active nematic particles we use the tensor order parameter $Q_{\alpha\beta}^{j}=\left(n_{\alpha}^{j}n_{\beta}^{j}-\delta_{\alpha\beta}/3\right)$, where the upper index corresponds to the particles and the lower indices represent the componenets $x,y,z$. The active force does not distinguish ’head from tail’ and it thus has the form:

[TABLE]

The torque reflects the fact that both parallel and anti-parallel configurations are favored. It has the form:

[TABLE]

The pair-interaction force $\mathbf{F}_{ij}$ and the holomonic constraint $g(\mathbf{q})$ are the same as in the active polar particles case. The simulation parameters for this case are $(J,k,\sigma,\gamma,\gamma_{a},v_{0})=(10,3,2,0.1,2.5,1.1)$ unless otherwise specified.

.4 Coarse-grained defect description

In the coarse-grained defects description by active polar particles Giomietal_PRL_2013 and Keberetal_Science_2014 the elastic energy between defects is $E\sim\log(q_{ij}^{g})$, where $q_{ij}^{g}$ is the geodesic distance between the defects. The pair-interaction force is therefore $\mathbf{F}_{ij}=\frac{k}{q_{ij}^{g}}\frac{\mathbf{q}_{i}-\mathbf{q}_{j}}{q_{ij}}$, which is no longer short-ranged. Defects align anti-parallel to each other and the restoring torque strength is Keberetal_Science_2014 $T_{i}=J\sum_{j\in U(i)}\cot(\frac{\theta_{ij}}{2})$, where $\theta_{ij}$ is the angle between $\mathbf{n}_{i}$ and $\mathbf{n}_{j}$. The vector form for the torque can be written in terms of the orientations as:

[TABLE]

Finally the defects are treated as self-propelled particles and the active force is $\mathbf{F}_{i}^{ac}=v_{0}\mathbf{n}_{i}$. The simulation parameters for this case are $(J,k,\gamma,\gamma_{a},v_{0})=(3,4,0.1,2.5,0.11)$ unless otherwise specified.

.5 Geometric properties

Besides a sphere we consider two classes of ellipsoidal surfaces: (i) spheroidal and (ii) non-spheroidal. These ellipsoids are characterized by their major axis $a$, $b$ and $c$ and have non-constant Gaussian curvature

[TABLE]

For spheroidal ellipsoids two of these values are equal. The algebraic description reads $g(\mathbf{q}_{i})=\frac{q_{i,1}^{2}}{a^{2}}+\frac{q_{i,2}^{2}}{b^{2}}+\frac{q_{i,3}^{2}}{c^{2}}-1=0$. An umbilic point is a point where the maximum and minimum curvatures coincide. At an umbilical point, the surface is ”locally spherical”. These points are found at

[TABLE]

In figure 7 we show three different ellipsoids, where umbilical points are highlighted and the color coding corresponds to the Gaussian curvature $K$.

.6 Numerical methods

Eqs. 1 have been numerically solved using RATTLE discretization HamDyn . The equation for the orientational dynamic eqs. 3 have been first solved unconstrained with the torque $\mathbf{T_{i}}$ projected on to the normal plane of the surface at point $\mathbf{r}_{i}$. Afterwards the orientation $\mathbf{n}_{i}$ has been projected on to the tangent plane of the surface at point $\mathbf{r}_{i}$ and the angular velocity $\bm{\omega}_{i}$ takes the direction of the normal to the surface at point $\mathbf{r}_{i}$.

We fix the number of particles $N=1000$ and the volume fraction $\phi\simeq 1$. The surface area for the sphere is equal to $A=4\pi R^{2}$, with $R=31.6$. The ellipsoid parameters $a,b,c$ have been chosen such that the surface area is equivalent to the surface area of the sphere and the aspect ratio is respected.

The nematic order parameter is defined as

[TABLE]

where the sum is over the nearest neighbors and $w_{ij}=q_{ij}^{-1}$.

Defects are calculated as the local center of mass for regions where the local order parameter $P_{i}$ is smaller then $0.45$ (some corrections were required for regions of high Gaussian curvature, due to strong distortion of the director field).

The simulation code is implemented in C++, using the GeographicLib library karney2011 for the calculation of the geodesic distances. However, for non-spheroidal ellipsoids the euclidian distance has been used. This approximation can be justified by the short-range interactions. Data have been analyzed using Python, Ovito ovito and Paraview.

Acknowledgements.

This work was funded by the European Union (ERDF) and the Free State of Saxony via the ESF project 100231947 (Young Investigators Group Computer Simulation for Materials Design - CoSiMa). We used computing resources provided by JSC within project HDR06.

Author contributions F.A. and C.K. derived the model and conducted the simulations, F.A. and A.V. analyzed the results. All authors contributed to a critical discussion of the results and participated in writing the manuscript, which A.V. finalized.