Curvature Potential Unveiled Topological Defect Attractors

: We consider the theoretical and positional assembling of topological defects (TDs) in effectively two-dimensional nematic liquid crystal ﬁlms. We use a phenomenological Helfrich– Landau–de Gennes-type mesoscopic model in which geometric shapes and nematic orientational order are expressed in terms of a curvature tensor ﬁeld and a nematic tensor order parameter ﬁeld. Extrinsic, intrinsic, and total curvature potentials are introduced using the parallel transport concept. These potentials reveal curvature seeded TD attractors. To test ground conﬁgurations, we used axially symmetric nematic ﬁlms exhibiting spherical topology.


Introduction
Topological defects (TDs) [1] refer to localized regions in an ordered manifold, where the local order is frustrated. They are present in all branches of physics because the sole condition for their existence is symmetry breaking [2]. Several key features of TDs could be extracted by studying two-dimensional (2D) systems, which could be efficiently mathematically and numerically approached [3]. Two-dimensional systems are, in particular, convenient as a means of analyzing the impact of geometry (curvature, topology) on the creation, stabilization, and positional assembling of TDs [4].
Two-dimensional systems hosting an orientational order field could exhibit point TDs [1]. Their key feature is characterized by a topological charge, which is a quantized and conserved quantity [5]. In 2D, this is equivalent to the winding number m [5], which elucidates the total reorientation of the respective order parameter field on encircling the defect counterclockwise. One commonly refers to TDs exhibiting m > 0 and m < 0 as defects and antidefects. In general, a nearby pair {defect, antidefect} = {m, −m} tends to annihilate into a local defectless state. Namely, TDs introduce local strong energetically expensive elastic distortions. For this reason, within defect cores, the ordering field is commonly melted or exhibits a qualitatively different structure with respect to bulk order [6]. TDs are, in general, rare in bulk equilibrium [2]. However, TDs might introduce new qualitative and effective properties into a system or can be exploited in several applications. Consequently, it is of interest to find efficient mechanisms, which enable the creation, stabilization, and manipulation of TDs.
In 2D systems, TDs could be efficiently controlled via the curvature of the ordered manifold. One commonly distinguishes between intrinsic and extrinsic curvature [7,8]. Note that several theoretical approaches express elastic distortions in terms of covariant derivatives [7,8]. These penalize deviations of the ordering field from local geodesics (the latter are analogs of straight lines in Euclidean space). In such approaches, the extrinsic curvature contributions, which are sensitive to how a 2D curved manifold is embedded in 3D, are discarded. Note that the energy associated with the extrinsic term, as described in [7,8], was actually considered in essentially the same form in studying biological membrane shapes [9][10][11][12][13][14] without TDs (the deviatoric term in Appendix C). A simple analysis [7] suggests that both intrinsic and extrinsic curvatures are in general present, featuring in elastic-free energies weighted by elastic constants of comparable magnitude. Furthermore, in several cases, the local impact of intrinsic and extrinsic contributions on TDs could be conflicting [7,8]. However, research on the mutual influence of curvature and TDs is relatively scarce.
An ideal system to study the physics of TDs is constituted by nematic liquid crystals (NLCs) [15,16]. In bulk equilibrium, NLCs are commonly described in terms of the nematic director field n, which exhibits head-to-tail symmetry, and tends to be uniformly aligned in bulk equilibrium. Due to their unique combination of liquid crystallinity, softness, optical transparency, and diversity, TDs could be relatively easily created, manipulated, and observed. Furthermore, thin and effectively 2D LC structures could be relatively easily experimentally realized, e.g., in nematic LC shells or biological membranes [17][18][19][20][21].
In this paper, we consider the role of combined intrinsic and extrinsic curvatures on structures in 2D nematic shells. We introduce the geometric potentials [21,22] that will reveal "hotspots", to which TDs are attracted. The structure of the paper is as follows. In the Methods section, we present our model, where we describe 2D configurations using both nematic and curvature fields, represented in tensorial form. Curvature potentials are introduced along with axially symmetric shapes, which we use to illustrate curvaturedriven phenomena. In the Results section, we present examples, which demonstrate the usefulness of geometric potentials to predict positional assembling of LCs on curved substrates. In the Appendices A-C, we describe some technical details.

Methods
We consider a two-dimensional curved manifold exhibiting a nematic in-plane orientational order. At the mesoscopic level, this is commonly described by the nematic director field [16], n, which is of unit length and exhibits the head-to-tail symmetry (i.e., |n| = 1, and the states ±n are physically equivalent). An equilibrium ground state of n is uniformly aligned along a single symmetry breaking direction in flat 2D manifolds. On the contrary, curved manifolds could host TDs due to topological reasons [5]. In the following, we present how one could predict the spatial positions of TDs for a given substrate geometry without solving the relevant Euler-Lagrange equilibrium equations. For illustration purposes, we henceforth limit the study to closed surfaces exhibiting spherical topology.

Intrinsic and Extrinsic Curvature Contributions
Orientational order on a curved manifold is in general controlled by two qualitatively different mechanisms: the intrinsic and extrinsic curvature contributions [7,8]. We first present their main roles using a minimal model, in which order is solely represented by n, which we parametrize by an angle θ as n = cos θe 1 + sin θe 2 . (1) Here, the unit vectors {e 1 , e 2 } determine the local principal curvature directions [23]. These define the outer local surface normal of a local surface patch: v = e 1 × e 2 .
We express the elastic free energy density in terms of a single positive constant K [16] f e = K 2 |∇ s n| 2 .
This contribution enforces spatially homogeneous orientational order. The operator ∇ s (•) = ∇(•)(e 1 ⊗ e 1 + e 2 ⊗ e 2 ) stands for the surface gradient [24], which deprives the "common" gradient ∇ of the leg along v. As {e 1 , e 2 } correspond to principal curvature directions, we can write [25] Here, κ g1 (κ g2 ) is the geodesic curvature and σ 1 (σ 2 ) the principal curvature along e 1 (e 2 ), respectively. It follows from Equation (3) that where stands for the curvature tensor and is referred to as the spin connection [23,24]. The invariants of C, trace and determinant, yield the mean curvature H and the Gaussian curvature G: For later use, we also define the deviatoric curvature measure [9] as The first and the second terms in the right-hand side of Equation (4) are typical representatives of intrinsic and extrinsic curvature contributions. Note that, in the single elastic constant approximation, both terms are present and are weighted by the same constant K.
The intrinsic term is the origin of ordering frustration in cases G = 0, which could be resolved by introducing TDs. On the other hand, the extrinsic term acts as a local ordering field, which is present if the principal curvatures are different. The key roles of these terms are summarized in Appendix A. There, we illustrate that a local surface patch possessing G > 0 (G < 0) attracts TDs with a positive (negative) topological charge. Furthermore, regions exhibiting E < 0 tend to expel TDs. Namely, in them a local ordering along a preferred orientation is enforced, which is incompatible with local nematic ordering of TDs.

Helfrich-Landau-Type Phenomenological Model
We next consider a more detailed model in which we focus on the impact of curvature on the position of TDs in nematic orientational order. We introduce the geometric potentials derived from the model free energy, which efficiently determine the location of TDs for a given manifold geometry. Positions of TDs reflect the interplay between intrinsic and extrinsic curvature contributions.
We describe the nematic order by the tensor order parameter field Q. In the local curvature principal frame (see Equation (5)), we use the following parametrization [24] Q = q 0 (e 1 ⊗ e 1 − e 2 ⊗ e 2 ) + q m (e 1 ⊗ e 2 + e 2 ⊗ e 1 ).
Here, q 0 and q m are variational variables. In its eigenframe, determined by the nematic director field n (see Equation (1)) it can be expressed as [24]: Here, {n, n ⊥ } are the eigenvectors of Q, corresponding to eigenvalues {λ, −λ}, where λ ∈ [0, 1/2] reveals the amplitude of nematic order in an infinitesimally small surface patch determined by the surface normal v = n × n ⊥ .
In the spirit of the classical Landau-type approach, we expand the free energy density f = f H + f c + f e in terms of invariants constructed using Q and C. The classical Helfrich curvature ( f H ), condensation ( f c ), and elastic ( f e ) contributions are expressed as [21,26] We introduced only the most essential terms to explain key mechanisms controlling positions of TDs. The classical Helfrich vesicle curvature contribution [26] f H is weighted by a positive bending modulus κ; it describes resistance to the manifold bending deformations. The quantities α 0 , β, T * in f c are positive phenomenological constants. The condensation term enforces nematic orientational order below the critical temperature T c . In a flat geometry, T c = T * and the condensed bulk equilibrium nematic order for T < T * is given by In the elastic free energy term, we included the simplest symmetry, which allowed for intrinsic and extrinsic curvature contributions, which are weighted by positive intrinsic (k i ) and extrinsic (k e ) elastic moduli.

Scaling and Dimensionless Free Energy Terms
In our illustrations we restrict to closed 2D manifolds hosting in-plane nematic order, where T < T c . We consider axially symmetric geometries of surface area A, exhibiting spherical topology.
For mathematical and numerical convenience, we scale the tensor order parameter with respect to the reference bulk equilibrium order λ 0 (i.e., within a flat uniformly aligned nematic film, see Equation (12)). The curvature tensor and spatial coordinates are scaled with respect to describing the radius of a spherically shaped manifold of surface area A.
In addition to the geometrically imposed length scale R, an important role is also played by the nematic order parameter correlation length ξ. It describes a typical distance on which a locally perturbed nematic amplitude is recovered. We estimate it by The resulting dimensionless free energy density ( f → f R 2 /k i ) reads as where µ = k e λ 0 k i . Note that the extrinsic term is weighted against the intrinsic term by a dimensionless coefficient µ ∝ λ −1 0 , which tends to diverge on approaching T c from below. In the numerical analysis, we used the parametrization of Q given by Equation (9) in terms of the scaled variational fields {q 0 /λ 0 , q m /λ 0 }. We consider axially symmetric 2D closed manifolds referred to as a 3D Cartesian system (x, y, z), defined by unit vectors {e x , e y , e z }. We represent the position vector r defining the 2D shapes by where ρ(s) and z(s) are the coordinates of the shape profile in the (ρ, z)-plane, ϕ ∈ [0, 2π] stands for the azimuthal angle and s represents the arc length of the profile curve. We calculate 2D shapes and the corresponding nematic structures by minimizing the total free energy for a given relative volume Here, V refers to the volume enclosed by an axially symmetric shape of total surface area A, and V 0 = 4πR 3 /3 is the volume of the sphere having the same surface area, where R is given by Equation (13). Technical details are given in Appendix B.

Curvature Potentials
We introduce geometric curvature potentials [22] which reveal attracting or repelling regions for TDs within a curved manifold. These potentials are calculated for a given geometry of the manifold.
To this purpose, we need to identify local ground states of Q for a given manifold geometry. In flat geometry, this corresponds to a spatially homogeneous alignment of n along a single symmetry breaking direction, and λ = λ 0 , yielding g i = g e = 0. In curved manifolds the local ground state might exhibit a finite elastic penalty, which we refer to as the fossil elastic energy (see also [27]).
To determine this, we request that the nematic director is locally parallel transported [24]. In this case, it exhibits locally minimal elastic distortions. A locally parallel transported unit vector e (p) (the superscript (p) indicates that a vector is parallel transported) obeys the equation [24] ∇ s e (p) = −v ⊗ Ce (p) .
The corresponding scaled parallel transported nematic order tensor is expressed as Note that the assumption that λ = λ 0 is here understood. Taking into account Equations (1), (18) and (19) we obtain The equations above introduce the intrinsic curvature potential w i and the extrinsic curvature potential w e . One sees that w i is independent of θ. Next, w e is directly proportional to the deviatoric free energy term which was originally introduced in studying biological membrane shapes [12,14,28,29] (see Appendix C). For a positive value of k e (µ > 0), the extrinsic curvature potential tends to align n along the principal direction exhibiting minimal absolute curvature, for which w (min) e = − σ 2 1 − σ 2 2 . This term acts similar to an ordering field, which is incompatible with a local nematic TD. Consequently, regions, where the extrinsic curvature is present, do not favor the presence of TDs.
An additionally important parameter affecting the location of TDs is the effective local temperature. Namely, the presence of TDs in 2D requires local melting of nematic order. In our model, λ = 0 at the center of TDs. However, this is energetically costly below T c . Namely, in bulk equilibrium, where T c = T * , the condensation of nematic order is energetically advantageous and f (eq) = f Here, the superscript (eq) denotes the bulk equilibrium condition. The linear size of the defect core is roughly given by the nematic correlation length. Therefore, the free energy penalty ∆F for introducing a point defect is roughly given by Furthermore, elastically distorted regions exhibit effectively different temperatures. To illustrate this, we limit ourselves to consider the intrinsic curvature term, which is quadratic in λ 2 . Neglecting spatial variations in λ, and taking into account only quadratic contribution in λ, we have that where Therefore, in elastically distorted regions, the critical temperature is effectively reduced, which locally softens the nematic order. Consequently, such regions attract cores of TDs. Note that both geometric potentials are, in general, present in a distorted region. Their common impact on the local nematic amplitude is measured by which we refer to as the total curvature potential. Here, w (min) e corresponds to the minimal value of w e upon varying θ. Our approximate analysis suggests that TDs will be attracted to regions where w t > 0 exhibits a maximum. That is, in these regions, the nematic order is softer due to an effectively lower local phase transition temperature.

Results
Of our interest here is the impact of intrinsic and extrinsic curvature on positions of TDs in nematic order within axially symmetric closed surfaces which are accessible in our model. It is already well known that regions exhibiting G > 0 (G < 0) attract TDs bearing positive (negative) topological charge if the extrinsic curvature is neglected. We aim to primarily predict positions of TDs for a given substrate geometry, where both curvature contributions are taken into account. Figure 1 shows the typical, qualitatively different, and axially symmetric structures considered in our study and their essential geometric properties, i.e., a spatially varying Gaussian curvature G (Equation (7b)) and a nonvanishing deviatoric curvature E (Equation (7c)), both expressed as functions of the arc length s along the generating profile. They are commonly dubbed stomatocytes, oblates, and prolates, and are depicted in black, red, and blue color, respectively. The key parameter determining their stability is the reduced volume v.  . The diagrams reveal a curvature-preferred local orientation of . Furthermore, regions with a large enough | | do not favor the presence of TDs. In Figure 3, we present the total curvature potential for these shapes for varying s and , which measures the relative importance of extrinsic and intrinsic curvature contributions. One sees that has a particularly strong impact on oblate shapes (Figure 3b). Typical equilibrium configurations (shapes and nematic textures) in the absence ( = 0) and for a relatively strong presence of extrinsic curvature contribution ( = 1) are plotted in Figures 4 and 5, respectively. In the first panels, we plot nematic configurations in the ( , ) plane. Centers of defects are revealed by regions where = 0 and we mark them with capital letters A, B. The linear size of defect cores, where the nematic order is essentially melted, is roughly given by . In the second panels, we plot shapes of closed nematic shells, where we also indicate the location of TDs. In the third panels, we plot = ( ) dependence of corresponding shapes. One sees that maxima of ( ) are indeed attractors of TDs. The diagrams reveal a curvature-preferred local orientation of n. Furthermore, regions with a large enough |w e | do not favor the presence of TDs. In Figure 3, we present the total curvature potential w t for these shapes for varying s and µ, which measures the relative importance of extrinsic and intrinsic curvature contributions. One sees that µ has a particularly strong impact on oblate shapes (Figure 3b). Typical equilibrium configurations (shapes and nematic textures) in the absence (µ = 0) and for a relatively strong presence of extrinsic curvature contribution (µ = 1) are plotted in Figures 4 and 5, respectively. In the first panels, we plot nematic configurations in the (s, ϕ) plane. Centers of defects are revealed by regions where λ = 0 and we mark them with capital letters A, B. The linear size of defect cores, where the nematic order is essentially melted, is roughly given by ξ.
In the second panels, we plot shapes of closed nematic shells, where we also indicate the location of TDs. In the third panels, we plot w t = w t (s) dependence of corresponding shapes. One sees that maxima of w t (s) are indeed attractors of TDs.        Below, we discuss in more detail the role of curvature fields in determining the location of TDs. Note that, according to the Gauss-Bonnet and Poincaré-Hopf theorems [23], the total topological charge of TDs equals = 2 for the spherical geometry (see Equation TDs. On the contrary, the stomatocytes possess a region where < 0, which attracts < 0. Furthermore, | | is large enough in this region to trigger the nucleation of additional defect-antidefect pairs = , = − aiming to establish topological neutrality in each surface patch characterized by different values of the average Gaussian curvature (see Appendix A). For = 0 the positions of defects are dominated by intrinsic curvature and are attracted to maxima of ( ). Note that, in the case of prolate shapes, the defects are not exactly placed at maxima of ( ), because these points attract two TDs of the same charge and are, therefore, mutually repealed. On the contrary, in oblate structures, the equatorial region is large enough to accommodate all four TDs. In stomatocytes, antidefects and defects are attracted to the respective maxima of w (s) in regions exhibiting G < 0 and > 0, respectively. This is in line with the ETCC mechanism [20] (see Appendix Below, we discuss in more detail the role of curvature fields in determining the location of TDs. Note that, according to the Gauss-Bonnet and Poincaré-Hopf theorems [23], the total topological charge of TDs equals m = 2 for the spherical geometry (see Equation (A2) in Appendix A). Furthermore, highly charged TDs are energetically unfavorable for structures characterized by R ξ 1. For this reason, TDs bearing elementary charges |m| = 1/2 are favored. Consequently, prolate and oblate structures possess four m = 1/2 TDs. On the contrary, the stomatocytes possess a region where G < 0, which attracts m < 0. Furthermore, |G| is large enough in this region to trigger the nucleation of additional defect-antidefect pairs m = 1 2 , m = − 1 2 aiming to establish topological neutrality in each surface patch characterized by different values of the average Gaussian curvature G (see Appendix A).
For µ = 0 the positions of defects are dominated by intrinsic curvature and are attracted to maxima of w t (s). Note that, in the case of prolate shapes, the defects are not exactly placed at maxima of w t (s), because these points attract two TDs of the same charge and are, therefore, mutually repealed. On the contrary, in oblate structures, the equatorial region is large enough to accommodate all four TDs. In stomatocytes, antidefects and defects are attracted to the respective maxima of w t (s) in regions exhibiting G < 0 and G > 0, respectively. This is in line with the ETCC mechanism [20] (see Appendix A) and also the proposed "attractiveness" of w t (s) maxima. Note that, in G > 0, the regions' maxima are relatively shallow. Consequently, locations of m = 1/2 TDs are dominated by their mutual repulsion.
For µ = 1 we observe only quantitative changes in prolate shapes and additional qualitative changes in the remaining competitive shapes. In the prolate configuration, a relatively strong extrinsic curvature field appears outside the poles of the shape. Consequently, the four TDs are pushed closer to the poles in comparison to the µ = 0 case. In this case, the extrinsic and intrinsic curvatures impose similar preferences to TDs. On the contrary, in oblate shapes the intrinsic curvature favors positioning TDs in the equatorial region, while the extrinsic one tends to expel them from there. One sees that they are located at local maxima in w t (s). Note that they are not located at the global maximum because, there, two mutually repelling TDs would be located at the same point. Finally, in the stomatocyte configuration shown, the extrinsic term is strong enough to prevent the formation of additional defect pairs. Note that the region at s/L s ∼ 0.6 does not possess TDs. However, the nematic order there is strongly suppressed due to relatively strong local elastic distortions.
Finally, in Figure 6 we plot the stability diagram of the competing structures on varying v in the presence and absence of the extrinsic contribution. One sees that the latter strongly favors oblate-type shapes. A) and also the proposed "attractiveness" of ( ) maxima. Note that, in > 0, the regions' maxima are relatively shallow. Consequently, locations of = 1/2 TDs are dominated by their mutual repulsion. For = 1 we observe only quantitative changes in prolate shapes and additional qualitative changes in the remaining competitive shapes. In the prolate configuration, a relatively strong extrinsic curvature field appears outside the poles of the shape. Consequently, the four TDs are pushed closer to the poles in comparison to the = 0 case. In this case, the extrinsic and intrinsic curvatures impose similar preferences to TDs. On the contrary, in oblate shapes the intrinsic curvature favors positioning TDs in the equatorial region, while the extrinsic one tends to expel them from there. One sees that they are located at local maxima in ( ). Note that they are not located at the global maximum because, there, two mutually repelling TDs would be located at the same point. Finally, in the stomatocyte configuration shown, the extrinsic term is strong enough to prevent the formation of additional defect pairs. Note that the region at / ~0.6 does not possess TDs. However, the nematic order there is strongly suppressed due to relatively strong local elastic distortions. Finally, in Figure 6 we plot the stability diagram of the competing structures on varying in the presence and absence of the extrinsic contribution. One sees that the latter strongly favors oblate-type shapes.

Conclusions
Of interest in this study was the identification of attractor sites for topological defects (TDs) on two-dimensional (2D) curved substrates hosting nematic orientational order. In general, the positions of TDs are influenced by both intrinsic and extrinsic curvature contributions. Their impact on TDs may be antagonistic with regard to their geometries. Furthermore, in most studies so far, the role of extrinsic curvature contributions has been (in Figure 6. Stability diagram of equilibrium closed membrane shapes as a function of relative volume v in the presence (k e = k i /2, solid line) and in the absence (k e = 0, dashed line) of the extrinsic contribution. Profiles of stable shapes presented in the diagram for different values of the reduced volume v and k i /к were calculated in the presence of the extrinsic term (k e = k i /2). Stomatocyte shapes are colored in black, oblate shapes (including discocyte shapes) in red and prolate shapes in blue. R/ξ = 7.

Conclusions
Of interest in this study was the identification of attractor sites for topological defects (TDs) on two-dimensional (2D) curved substrates hosting nematic orientational order. In general, the positions of TDs are influenced by both intrinsic and extrinsic curvature contributions. Their impact on TDs may be antagonistic with regard to their geometries. Furthermore, in most studies so far, the role of extrinsic curvature contributions has been (in several cases unjustifiably) neglected and the combined effects of both curvature contributions are so far scarcely explored. In our study, we intended to contribute to a general understanding of their mutual influence on TDs in orientational order.
In our theoretical analysis, we described the local geometry by the curvature tensor and the nematic order by the nematic tensor order parameter. Positions of centers of points of TDs in our study were fingerprinted by points where the nematic order parameter amplitude was locally melted. To predetermine positions attracting TDs, we introduced intrinsic w i and extrinsic w e geometric potentials by applying the classical notion of parallel transport. This enabled us to determine local ground states for a given substrate geometry.
The intrinsic potential is independent of the local nematic director orientation θ, while the extrinsic term depends on θ. The regions characterized by |w e | 1 tend to expel TDs because, in them, geometrically enforced ordering is imposed, which is incompatible with locally relatively variable nematic structures of TDs. The minimization of w e with respect to θ yields its ground state value w (min) e . Our analysis reveals that maxima of the total geometric potential w t = w i + µ w (min) e are attracting TDs. Here, µ = k e /(k i λ 0 ) measures the relative importance of w i and w (min) e . Note that µ ∝ 1/λ 0 and, in general, one expects comparable values of k e and k i . Therefore, just below the nematic-isotropic phase transition the extrinsic contribution is expected to play a dominant role. We tested the predicting power of geometric potentials by studying curvature-imposed positional assembly of TDs in axially symmetric nematic shells, exhibiting spherical symmetry. Upon varying the relative volume v, the shells could exhibit three qualitatively different geometries, and in each TDs were attracted to maxima of w t , confirming our expectation.
Our study effectively illustrates that curved substrates could be exploited for efficient assembling of TDs to predetermined positions. This is advantageous for several applications. For example, the cores of fixed TDs could be exploited to manipulate the polarization of light [30]. Furthermore, TDs could efficiently trap appropriate nanoparticles [31,32]. Note that TDs in liquid crystals could be relatively easily manipulated [33][34][35][36][37][38][39][40][41], which offers an indirect path to manipulate assemblies of trapped NPs. Finally, understanding the curvature-enabled stabilization mechanisms of TDs and their assemblies might shed light on still unresolved problems in fundamental physics [42]. For example, it seems that fields represent fundamental natural entities, and topologically protected localized field distortions, such as TDs, might be related to "particles", if nature is viewed from this alternative perspective [43].

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
We first illustrate that the intrinsic curvature term can generate TDs on surface patches characterized by G = 0. Namely, its contribution K|∇ s θ + A| 2 in Equation (4) is removed if ∇ s θ = −A. Applying ∇× operation to this equation yields a frustrating condition which could be resolved by the local melting of order (i.e., by introducing a topological defect). Furthermore, according to the Gauss-Bonnet and Poincaré-Hopf theorems the integrated Gaussian curvature over a closed surface hosting in-plane order is quantized: Here, m stands for the total topological charge within the closed surface of genus g. The latter quantity counts the number of holes within the surface. For example, in spherical topology, g = 0 and m = 2. Equation (A2) suggests that one could imagine that a positive (negative) Gaussian curvature effectively acts as a smeared negative (positive) topological charge, which is neutralized by "real" discrete topological charge, carried by TDs. In this perspective one introduces the smeared curvature charge to a surface elementary patch d 2 r. For a closed surface, it strictly holds that m + m G = 0, so that the total system is topologically neutral (i.e., the total topological charge of the whole system equals zero). Hence, Equation (A3) suggests that a local region exhibiting [33] G > 0 (G < 0) attracts TDs bearing m > 0 (m < 0). This phenomenon is embodied in the Effective Topological Charge Cancellation (ETCC) mechanism [20,33], which refers to surfaces exhibiting spatially varying G. It claims that each surface patch ∆ζ, to which one can assign a characteristic average Gaussian curvature G[∆ζ] = 1 ∆ζ ∆ζ Gd 2 r, tends to be topologically neutral. Therefore, the total topological charge m[∆ζ] of TDs within it tends to compensate the corresponding total smeared curvature charge m G [∆ζ] = ∆ζ dm G .
For K > 0, this term is minimized for n aligned along the principal direction exhibiting the minimal curvature.

Appendix B
Here, we present technical details for calculating the equilibrium of 2D shapes and the corresponding nematic structures. Functions describing the coordinates of the shape profile ρ(s) and z(s) (see Equation (17)) are calculated as [20,[44][45][46]: where the function θ(s) represents the angle of the tangent to the profile curve with the plane that is perpendicular to the axis of rotation e z . For closed and smooth surfaces, we have to apply the following boundary conditions: θ(0) = 0, θ(L s ) = π, ρ(0) = ρ(L s ) = 0.
Where L s stands for the length of the profile curve. In our simulations, the function θ(s) is approximated by the Fourier series [20,[44][45][46]: where N is the number of Fourier modes, a i are the Fourier amplitudes, and θ 0 = θ(L s ) = π is the angle at the north pole of the surface. The local principal curvatures σ 1 and σ 2 , which appear in different energy contributions, are calculated as dθ(s) ds and sin θ(s) ρ(s) , respectively [20]. Our goal is to determine equilibrium closed shell shapes and the nematic ordering profiles on these shells. Equilibrium shell shapes are determined by the numerical minimization of a function of many variables [44][45][46]. These variables include the Fourier amplitudes a i and the shape profile length L s (see Equation (A6)). During the minimization procedure, the reduced volume v is kept constant. Next, we determine the equilibrium nematic configuration on previously calculated shape. The nematic order tensor is given by Equation (9), where q 0 and q m are variational parameters. Equilibrium configurations for q 0 and q m are calculated on a fixed shape by using the standard Monte Carlo method. The closed surface in the coordinates (ϕ, s) is represented as a network of 101 × 101 points. After that, the equilibrium shape is adjusted according to the current nematic texture. We repeat this process several times to obtain the equilibrium shell shape and the equilibrium nematic configuration.

Appendix C
Here, we demonstrate the similarity between the extrinsic term used in our study and the deviatoric term used in the study of biological membranes. Phospholipid bilayers forming a closed surface of biological membranes are assumed to be covered with anisotropic components oriented at different directions [25,26,47,48]. The orientation of the inclusion can be described by θ, which in this case represents the rotation of the principal directions of the inclusion's intrinsic shape relative to the membrane principal directions. The energy of a single inclusion is proportional to the mismatch between the intrinsic shape of the inclusion and the local membrane shape [9,11,13,26,47,48]: where K 1 and K 2 are constants, H is the local mean curvature of the membrane given by Equation (7a) and D = (σ 1 − σ 2 )/2 is the curvature deviator of the membrane. Inclusion's intrinsic shape is described by its intrinsic mean curvature H m and its intrinsic curvature deviator D m [25,26,47,48]. If D m = 0 the inclusion is isotropic while if D m = 0 it is anisotropic. In this formulation, membrane components are two-dimensional (described by two parameters). For minimal elastic distortions, the energy associated with the extrinsic term used in this study can be written as (see Equations (1), (18), (19) and (20b)): where k e is the extrinsic elastic modulus. We notice that the deviatoric term (Equation (A7)) and the extrinsic term (Equation (A8)) have the same dependence on the orientation of inclusions proportional to cos(2θ). The main difference is that the deviatoric term takes into account also the anisotropy of the inclusion's intrinsic shape described by D m . If we focus only on the orientation dependent part of the energy, we can write where γ is a constant, including different constants and curvatures, which can be determined by comparing Equation (A9) with Equations (A7) and (A8). With the aid of statistical mechanics, we can calculate the average orientation of a single inclusion [47]: where k B is the Boltzmann constant and T is the temperature. The result of the integration includes the quotient of the modified Bessel functions I 0 and I 1 [47]: The average orientation of the inclusion as a function of γ is presented in Figure A1 [47].
where is the Boltzmann constant and is the temperature. The result of the integration includes the quotient of the modified Bessel functions and [47]: The average orientation of the inclusion as a function of is presented in Figure A1 [47]. It can be seen in Figure A1 that the orientational ordering increases with increasing . In case of Equation (A7), is proportional to the inclusion's intrinsic curvature deviator and the curvature deviator of the membrane . The effect of orientational ordering is in this case strong in highly curved anisotropic regions containing anisotropic membrane components. In case of Equation (A8), membrane inclusions are modelled as straight rods and is proportional to the mean membrane curvature and the curvature deviator of the membrane . Strong orientational ordering is, therefore, present in highly curved anisotropic membrane regions except on saddle-like surfaces for which = 0 and ≠ 0.
For one-dimensional curved membrane components, we can use a model proposed in [49][50][51][52]. The elastic energy of a flexible rod-like membrane component (e.g., curved protein) is given by: where is a constant, is the intrinsic curvature of the crescent-shaped membrane domain and is the local membrane curvature seen by the attached domain, which can be expressed by Euler relation as = + cos (2 ). Here, is the angle between the orientation of the one-dimensional membrane domain and the principal curvature . Using this model, it was shown that curved BAR protein domains can induce the growth of tubular protrusions [49,52]. Without the application of external force, BAR domains were found to be oriented perpendicular to the protrusion in order to minimize the elastic energy given by Equation (A12) [52]. This theoretical prediction was recently confirmed in [53,54]. Equation (A12) can also be written as: We can observe that some parts of Equation (A13) are the same as in the extrinsic term given by Equations (A7) and (A8). It can be seen in Figure A1 that the orientational ordering increases with increasing γ. In case of Equation (A7), γ is proportional to the inclusion's intrinsic curvature deviator D m and the curvature deviator of the membrane D. The effect of orientational ordering is in this case strong in highly curved anisotropic regions containing anisotropic membrane components. In case of Equation (A8), membrane inclusions are modelled as straight rods and γ is proportional to the mean membrane curvature H and the curvature deviator of the membrane D. Strong orientational ordering is, therefore, present in highly curved anisotropic membrane regions except on saddle-like surfaces for which H = 0 and D = 0.
For one-dimensional curved membrane components, we can use a model proposed in [49][50][51][52]. The elastic energy of a flexible rod-like membrane component (e.g., curved protein) is given by: where κ is a constant, C p is the intrinsic curvature of the crescent-shaped membrane domain and C is the local membrane curvature seen by the attached domain, which can be expressed by Euler relation as C = H + D cos(2θ). Here, θ is the angle between the orientation of the one-dimensional membrane domain and the principal curvature σ 1 . Using this model, it was shown that curved BAR protein domains can induce the growth of tubular protrusions [49,52]. Without the application of external force, BAR domains were found to be oriented perpendicular to the protrusion in order to minimize the elastic energy given by Equation (A12) [52]. This theoretical prediction was recently confirmed in [53,54]. Equation (A12) can also be written as: W p (θ) = κ H 2 + 2HD cos(2θ) − 2HC p − 2DC p cos(2θ) + D 2 cos 2 (2θ) + C p 2 . (A13)