Annihilation of Highly-Charged Topological Defects

: We studied numerically external stimuli enforced annihilation of a pair of daughter nematic topological defect (TD) assemblies bearing a relatively strong topological charge | m | = 3/2. A Landau-de Gennes phenomenological approach in terms of tensor nematic order parameter was used in an e ﬀ ectively two-dimensional Cartesian coordinate system, where spatial variations along the z-axis were neglected. A pair of { m = 3/2, m = − 3/2} was enforced by an appropriate surface anchoring ﬁeld, mimicking an experimental sample realization using the atomic force microscope (AFM) scribing method. Furthermore, defects were conﬁned within a rectangular boundary that imposes strong tangential anchoring. This setup enabled complex and counter-intuitive annihilation processes on varying relevant parameters. We present two qualitatively di ﬀ erent annihilation paths, where we either gradually reduced the relative surface anchoring ﬁeld importance or increased an external in-plane spatially homogeneous electric ﬁeld E . The creation and depinning of additional defect pairs (cid:110) 12 , − 12 (cid:111) mediated the annihilation in such a geometry. Furthermore, we illustrate the absorption of TDs by sharp edges of the conﬁning boundary, accompanied by m = ± 1/4 ↔ ∓ 1/4 winding reversal of edge singularities, and also E -driven zero-dimensional to one-dimensional defect core transformation.


Introduction
Topological defects (TDs) [1] correspond to localized regions within ordered media where the relevant order parameter field is frustrated and topologically protected. In most cases, at the defect origin, the continuum field, which describes defectless equilibrium configurations, is not uniquely defined [2]. In general, TDs appear as a consequence of the universal concept of broken symmetry and are of interest in all branches of physics [1,3]. The key topological feature of a topological defect is fingerprinted in its topological charge [1,2], which is a conserved quantity. Corresponding topological charge conservation laws regulate transformations (decomposing, merging, depinning) of defect configurations.
Owing to their topological origin, the physics of TDs exhibits several universal features. Consequently, it is of interest to find physical systems in which TDs are relatively easily experimentally formed, manipulated, and observed, which allows detailed and controlled analysis of the behavior of interest. For example, in the last decades, several studies analyzed quench driven coarsening of TDs in different condensed matter systems to gain information on structural details in the early universe [3][4][5].
Nematic liquid crystals (NLC) [6] offer an ideal testbed to study TDs [7][8][9][10] due to their unique combination of liquid behavior, softness, experimentally-tractable characteristic scales of time and spatial responses, and optical anisotropy. Furthermore, NLC structures possessing TDs could be exploited in various (in particular electro-optical) applications [11,12]. The uniaxial nematic structure is commonly described by the mesoscopic nematic director field , which points along the local uniaxial NLC order, and the states ± are physically equivalent (the so-called head-to-tail invariance) [7]. In bulk, nematic equilibrium is homogeneously aligned along a symmetry breaking direction. Several studies have been performed on thin LC films, in which relatively simple quasi-two-dimensional (2D) behavior [13][14][15][16][17] is observed. Despite the simplicity, such systems display rich and complex behavior, in particular with respect to topological defect transformations. In 2D NLCs, only point defects exist. The role of topological charge is played by the winding number m [7] (also referred to as the Frank index). It determines the total rotation of on encircling the defect core counterclockwise. Due to the head-to-tail invariance, m can possess half-integer and whole integer numbers. Some typical TDs are shown in Figure 1. TDs bearing positive and negative m are commonly referred to as defects and antidefects, respectively, and pairs { , − } in general tend to annihilate into a defectless configuration. In general, high values of m are rarely realized. Namely, the free energy costs of a defect scales with m 2 [16]. Consequently, if one enforces locally a defect bearing | | > 1/2, it tends to decompose [16][17][18] to elementary defects bearing | | = 1/2, where the charge conservation law is obeyed. Note that the core structure of single | | = 1/2 disclination is well known. In terms of the tensor nematic order parameter it was originally determined by Schopohl and Sluckin [13]. Typical structural changes by traversing the defect core are shown in Figure 2. The orientational frustration is resolved via the order reconstruction (OR) mechanism [19]. In general, this mechanism is realized in cases where a relatively large orientational mismatch in nematic order is imposed on a scale comparable [19][20][21][22][23] to the nematic biaxial order parameter relaxation length . The latter is relatively weakly temperature-dependent and is of order a few tens of nanometer [23]. The presence of the order reconstruction (OR) mechanism within the core was confirmed experimentally [24] in lyotropic In general, high values of m are rarely realized. Namely, the free energy costs of a defect scales with m 2 [16]. Consequently, if one enforces locally a defect bearing |m| > 1/2, it tends to decompose [16][17][18] to elementary defects bearing |m| = 1/2, where the charge conservation law is obeyed. Note that the core structure of single |m| = 1/2 disclination is well known. In terms of the tensor nematic order parameter it was originally determined by Schopohl and Sluckin [13]. Typical structural changes by traversing the defect core are shown in Figure 2. The orientational frustration is resolved via the order reconstruction (OR) mechanism [19]. In general, this mechanism is realized in cases where a relatively large orientational mismatch in nematic order is imposed on a scale comparable [19][20][21][22][23] to the nematic biaxial order parameter relaxation length ξ b . The latter is relatively weakly temperature-dependent and is of order a few tens of nanometer [23]. The presence of the order reconstruction (OR) mechanism within the core was confirmed experimentally [24] in lyotropic chromonic LCs. Namely, such LCs exhibit a micron-scale nematic biaxial order parameter correlation length ξ b that provides experimental insight into the core structure, whose characteristic linear size is comparable to ξ b . Crystals 2020, 10, x FOR PEER REVIEW  3 of 13 chromonic LCs. Namely, such LCs exhibit a micron-scale nematic biaxial order parameter correlation length that provides experimental insight into the core structure, whose characteristic linear size is comparable to . Figure 2. Core structure of a typical m=1/2 disclination. (a) The biaxial spatial variation, measured via the biaxiality parameter (see Equation (3)). The core structure is characterized by a rim, where = 1. The rim is circular in the approximation of equal elastic constants. The center of the core exhibits negative uniaxiality. (b) Characteristic changes in the ( , ) phase space on crossing the center of disclination. (c) Schematic variation of the nematic mesoscopic local profile within the core. Here numbers 1 to 5 mark points representing mesoscopic changes on traversing the defect core along a chosen direction.
In this paper, we focus on the annihilation of surface enforced relatively highly charged {defect, antidefect} pair bearing the winding number | | = 3 2 ⁄ in an effectively 2D nematic cell confinement.
We enforce the annihilation via different routes: i) by changing effective elastic LC properties, or ii) by increasing an external in-plane electric field strength. We demonstrate that annihilation proceeds via two qualitatively different channels, whose realization is rather counterintuitive.

Methods
We use the Landau-de Gennes mesoscopic approach in terms of the traceless and symmetric tensor nematic order parameter [6]. In its eigen frame, it is expressed as [20] = ⨂ , where and stand for eigenvalues and eigenvectors, respectively. Uniaxial states are commonly described by [6] = ⨂ − .
The nematic uniaxial order parameter ∈ [−1/2,1] quantifies the extent of fluctuations about the local nematic director , and > 0 ( < 0) reflects "prolate" ("oblate") shape in the ellipsoidal presentation of the LC mesoscopic order. We henceforth refer to such configurations as positive uniaxial and negative uniaxial nematics.  (3)). The core structure is characterized by a rim, where β 2 = 1. The rim is circular in the approximation of equal elastic constants. The center of the core exhibits negative uniaxiality. (b) Characteristic changes in the (s, ψ) phase space on crossing the center of disclination. (c) Schematic variation of the nematic mesoscopic local profile within the core. Here numbers 1 to 5 mark points representing mesoscopic changes on traversing the defect core along a chosen direction.
In this paper, we focus on the annihilation of surface enforced relatively highly charged {defect, antidefect} pair bearing the winding number |m| = 3/2 in an effectively 2D nematic cell confinement. We enforce the annihilation via different routes: (i) by changing effective elastic LC properties, or (ii) by increasing an external in-plane electric field strength. We demonstrate that annihilation proceeds via two qualitatively different channels, whose realization is rather counterintuitive.

Methods
We use the Landau-de Gennes mesoscopic approach in terms of the traceless and symmetric tensor nematic order parameter Q [6]. In its eigen frame, it is expressed as [20] where λ i andê i stand for Q eigenvalues and eigenvectors, respectively. Uniaxial states are commonly described by [6] The nematic uniaxial order parameter S ∈ [−1/2, 1] quantifies the extent of fluctuations about the local nematic directorn, and S > 0 (S < 0) reflects "prolate" ("oblate") shape in the ellipsoidal Crystals 2020, 10, 673 4 of 13 presentation of the LC mesoscopic order. We henceforth refer to such configurations as positive uniaxial and negative uniaxial nematics.
Biaxial states could be assessed if elastic distortions are present. The degree of biaxiality is well measured by the biaxiality parameter [25] The limiting values β 2 = 0 and β 2 = 1 correspond to a uniaxial nematic state and a nematic order exhibiting the maximal biaxiality, respectively.
The eigenvalues can be parametrized as where Note that the eigenvalues {λ 1 , λ 2 , λ 3 } can be expressed using two independent variables (i.e., {s,ψ}) due to the constraint 3 i=1 λ 2 i = 0. The (s,ψ) order parameter space is shown in Figure 2b. The isotropic order is reflected by s = 0. Nematic configurations corresponding to ψ = 0, ψ = 2π 3 , and ψ = − 2π 3 exhibit positive uniaxial states alongê 1 ,ê 2 , andê 3 , respectively. States ψ = π, ψ = − π 3 , and ψ = π 3 correspond to the negative uniaxial states alongê 1 ,ê 2 , andê 3 . The remaining states are biaxial. In this parametrization it holds β 2 = sin 2 (3ψ). Any continuous change in nematic order is fingerprinted in a continuous change in the order parameter space. In our study we will consider nematic structures in the Cartesian coordinate system (x,y,z), whose coordinate unit vectors are given by the triad (ê x ,ê y ,ê y ). Of interest are quasi-two-dimensional structures, where the Q eigenvector along the z-axis remains fixed (i.e.,ê 3 =ê z ). In this case, in addition to {s,ψ}, only one additional parameter (i.e. an angle θ) is needed to determine Q, describing the rotation of the Q eigenframe with respect to the Cartesian frame: Presentation of Q in terms of {s,ψ, θ} is useful to illustrate or guess different possible Q configuration variations along a given path in real space connecting two sites imposing different nematic order. An example is illustrated in Figure 2, where one crosses the core of |m| = 1/2 defect.
However, this Q representation does not provide a one-to-one mapping of each component (is not injective) [26]. For this reason, we used in simulations parametrization given by in terms of variational fields {q 1 , q 2 , q 3 }. It holds λ 1 = q 3 + q 2 1 + q 2 2 , λ 2 = q 3 − q 2 1 + q 2 2 , and λ 3 = −2q 3 . Note that the exchange of eigenvalues [13] between λ 1 and λ 2 corresponds to the condition q 2 1 + q 2 2 = 0. It was originally introduced to describe the elastic deformation within the core structure of |m| = 1/2 wedge defects [13]. The basic mechanism enabling this transformation is commonly referred to as order reconstruction [12,19].

Free Energy
We express the free energy F of a confined LC as a sum of nematic condensation ( f c ), elastic ( f e ), external electric field f f , and surface ( f (i) s ) free energy density contributions: Here d 3 r and d 2 r stand for the volume and area measures, and the superscript (i) refers to the i-th confining substrate. The free energy densities are expressed using the standard Landau-de Gennes expansion [6] in Q. We take into account only the most essential symmetry allowed contributions to illustrate phenomena of interest [6,20,27]: Here a 0 , b, c are positive material-dependent constants, T * stands for the isotropic supercooling temperature, L is a positive representative nematic elastic constant (i.e., we use the single elastic constant approximation), E stands for an external electric field, ∆ε is the dielectric anisotropy (we limit to materials with ∆ε > 0), w (i) > 0 is the anchoring strength coefficient at the i-th confining surface, which enforces nematic order described by Q s [27].

Geometry of the Problem
We study nematic configurations inside a three-dimensional well with square cross section of length R and thickness h<<R. The geometry of the problem is depicted in Figure 3, where we use Cartesian coordinate system (x,y,z). At the lateral walls ( Figure 3a) we enforce strong uniaxial tangential anchoring conditions: At the bottom (command) plate we enforce the nematic pattern consisting of a pair {defect, antidefect} of strength {m = m 1 = −3/2, m = m 2 = 3/2} placed at (x 1 , y 1 ) and (x 2 , y 2 ): Typically, we set (x 1 = 0, y 1 = −0.25 R) and (x 2 = 0, y 2 = 0.25 R). The opposite plate enforces degenerate tangential anchoring.
This setup mimics samples that could be realized experimentally using, for instance, the atomic force microscope (AFM) scribing method [28]. In a typical experimental setup, a nematic LC is confined within a thin plane-parallel cell, where at least one (command) surface imposes anchoring conditions inscribed via an AFM stylus [28], while the opposing surface imposes a degenerate tangential anchoring.   In simulations, we assume that the cell is thin enough so that the nematic structure is effectively two dimensional, exhibiting only variations in the (x,y) plane. Consequently, we use in numerical simulation the parametrization for the nematic tensor order parameter [22] given by Equation (7), where q 1 = q 1 (x, y), q 2 = q 2 (x, y), and q 3 = q 3 (x, y) are the variational parameters. In this parametrisationê z is always an eigenvector of Q.
For presentational purposes, we introduce the dimensionless temperature t = (T − T * )/(T * * − T * ) and material characteristic lengths [20]: The quantity ξ f stands for the external field nematic extrapolation length, which we express at the superheating temperature T = T * * , ξ b is the biaxial correlation length, d e is the surface extrapolation length, w stands for the anchoring strength imposed by the command surface, and S * * = 3b/16c stands for the superheating bulk value of S. We henceforth utilize the relative anchoring strength of the command surface plate and of the external field E = Eê y in terms of dimensionless ratios, respectively.
We minimize the free energy of the system with respect to the variational parameters q 1 , q 2 , and q 3 . The resulting Euler-Lagrange equations are solved numerically using the standard over-relaxation method. Technical and numerical details are given in [29].

Results
Of interest are different paths in which annihilation of a master surface anchoring enforced pair {m = 3/2, m = −3/2} annihilates in the geometry and boundary conditions described in the Methods section. The initial nematic configuration is shown in Figure 4a (biaxiality pattern) and Figure 5a (director field pattern). Each of the command surface imposed |m| = 3/2 defects at (x 1 = 0, y 1 = −0.25 R) and (x 2 = 0, y 2 = 0.25 R) decay in three elementary |m| = 1/2 daughter defects. Due to boundary conditions, two additional defects m = 1/2 and m = −1/2 at the top and the bottom boundaries of Figure 4a, respectively, are introduced. In addition, at the lateral boundaries, relatively (non-singular) strongly elastically deformed regions exist, which are fingerprinted by relatively strong local biaxiality. Furthermore, sharp corners of their structure enforce |m| = 1/4 type surface distortions [30].
For latter convenience, we label the three positively charged daughter defects (emerging from m = 3/2) by P1, P2, P3, and the corresponding negatively charged daughter defects (emerging from m = −3/2) by N1, N2, N3 as shown in Figure 4a. Boundary condition generated defects are labeled by P4 (upper boundary) and N4 (bottom boundary). Furthermore, at some stages at lateral boundaries, elastic distortions are strong enough to trigger nucleation of creations of pairs {defect, antidefect}, see Figure 4c. We label these defect pairs by (N5,P5) and (N6,P6).   First, we enforce the annihilation by gradually decreasing the relative importance of the surface anchoring field. In practice, this could be achieved by decreasing temperature towards the second-order nematic-SmA phase transition temperature T NA , which is accompanied by divergence [6] in the nematic twist and bend elastic constants. The relative importance of elastic and surface penalties of a confined NLC is measured by the ratio RW/K, where R stands for the characteristic confinement length, and K ∼ LS 2 is the representative Frank nematic elastic constant, and W measures the dimensional anchoring strength. In our simulations this ratio is measured by µ w . Therefore, on approaching T NA , the ratio µ w would decrease. The limit µ w 1 corresponds to the strong anchoring regime, while µ w < 1 corresponds to weak anchoring. The initial configuration, shown in Figure 4a, is evaluated for µ w = 1, which is strong enough to prevent the annihilation of {3/2,-3/2} pair; however, weak enough to allow decomposition of |m| = 3/2 singularities. Such conditions are realized in typical samples prepared using AFM scribing method [17]. For latter convenience, we label the three positively charged daughter defects (emerging from = 3/2) by P1, P2, P3, and the corresponding negatively charged daughter defects (emerging from = −3/2) by N1, N2, N3 as shown in Figure 4a. Boundary condition generated defects are labeled by P4 (upper boundary) and N4 (bottom boundary). Furthermore, at some stages at lateral boundaries, elastic distortions are strong enough to trigger nucleation of creations of pairs {defect, antidefect}, see Figure 4c. We label these defect pairs by (N5,P5) and (N6,P6).
First, we enforce the annihilation by gradually decreasing the relative importance of the surface anchoring field. In practice, this could be achieved by decreasing temperature towards the secondorder nematic-SmA phase transition temperature , which is accompanied by divergence [6] in the nematic twist and bend elastic constants. The relative importance of elastic and surface penalties of a confined NLC is measured by the ratio ⁄ , where R stands for the characteristic confinement length, and ~ is the representative Frank nematic elastic constant, and W measures the dimensional anchoring strength. In our simulations this ratio is measured by . Therefore, on approaching , the ratio would decrease. The limit ≫ 1 corresponds to the strong anchoring regime, while < 1 corresponds to weak anchoring. The initial configuration, shown in Figure 4a, is evaluated for = 1, which is strong enough to prevent the annihilation of {3/2,-3/2} pair; however, weak enough to allow decomposition of | | = 3 2 ⁄ singularities. Such conditions are realized in typical samples prepared using AFM scribing method [17]. On decreasing , in the first stage, the closest oppositely charged daughter defects (N1,P1) begin to approach (see Supplementary movie Annihilation1) until they mutually annihilate (See On decreasing µ w , in the first stage, the closest oppositely charged daughter defects (N1,P1) begin to approach (see Supplementary movie Annihilation 1) until they mutually annihilate (See Figure 4b).
Then the defect pairs (N2,P4), (P2,N4), and afterward (P3,N3) annihilate. During these processes, rearrangement of defects imposes sufficiently strong elastic distortions at the lateral boundaries to trigger the formation and depinning of two pairs of defects (P5,N5) and (P6,N6), as shown in Figure 4c. Finally, these defects are adsorbed by sharp edges of the confining boundary. Consequently, the winding number, characterizing strong distortions in corners, revert their winding number. This is evidently shown in Figure 5c,d, where one observes ±1/4 ↔ ∓1/4 winding number transformations of the director field profiles within the edges.
Next, we trigger annihilation by applying and gradually increasing a spatially homogeneous external electric field E = Eê y . The initial configuration for E = 0 is shown in Figure 6a, where we set (x 1 = 0, y 1 = −0.3 R) and (x 2 = 0, y 2 = 0.3 R). On increasing E (see Supplementary movie Annihilation 2), in the first stage, the inner (N1,P1) defect approach and finally annihilate (Figure 6c). The remaining defects annihilate in a qualitatively different way compared to the process shown in Figure 4, where lateral defect pairs (N5,P5) and (N6,P6) are strongly involved. Namely, the oppositely charged daughter and lateral defects approach each other and finally annihilate (Figure 6d). Furthermore, initially localized |m| = 1/2 structures of P4 and P5, clearly visible in Figure 6a-e, transform to nonlocalized order reconstruction planes. This defect core "explosion" was studied in detail in [31]. Here we summarize the key stages of the process. In this E-driven transformation, the initially spherically shaped biaxial core of a defect begins to grow as it is pushed towards a limiting wall (in our setting y = R/2 boundary for P4, and y = −R/2 boundary for M4) on increasing E. For E > 0, the ellipsoidal core-shape becomes progressively elongated along the x-axis. At some critical value of E, the core diverges (i.e., its longest linear size equals R), and a planar order reconstruction plane is formed. In this case order reconstruction takes form near the bottom and top confining plate. The corresponding nematic director profile, reflecting the Q eigen frame orientation, is shown Figure 7b. Note that on crossing the OR plane the LC molecules experience a mesoscopic shape transformation, similar to the one shown in Figure 2c along the direction 1-5. The OR plane is formed at µ f ∼ 3.5, which for typical LCs [6] corresponds to E ∼ 10 8 V/m.
Consequently, the winding number, characterizing strong distortions in corners, revert their winding number. This is evidently shown in Figure 5c,d, where one observes ±1/4 ↔ ∓1/4 winding number transformations of the director field profiles within the edges.
Next, we trigger annihilation by applying and gradually increasing a spatially homogeneous external electric field = . The initial configuration for E=0 is shown in Figure 6a, where we set ( = 0, = −0.3 ) and ( = 0, = 0.3 ). On increasing (see Supplementary movie Annihilation 2), in the first stage, the inner (N1,P1) defect approach and finally annihilate (Figure 6c). The remaining defects annihilate in a qualitatively different way compared to the process shown in Figure 4, where lateral defect pairs (N5,P5) and (N6,P6) are strongly involved. Namely, the oppositely charged daughter and lateral defects approach each other and finally annihilate (Figure 6d). Furthermore, initially localized | | = 1 2 ⁄ structures of P4 and P5, clearly visible in Figure 6a-e, transform to nonlocalized order reconstruction planes. This defect core "explosion" was studied in detail in [31]. Here we summarize the key stages of the process. In this E-driven transformation, the initially spherically shaped biaxial core of a defect begins to grow as it is pushed towards a limiting wall (in our setting y=R/2 boundary for P4, and y=-R/2 boundary for M4) on increasing E. For E>0, the ellipsoidal core-shape becomes progressively elongated along the x-axis. At some critical value of E, the core diverges (i.e., its longest linear size equals R), and a planar order reconstruction plane is formed. In this case order reconstruction takes form near the bottom and top confining plate. The corresponding nematic director profile, reflecting the Q eigen frame orientation, is shown Figure 7b. Note that on crossing the OR plane the LC molecules experience a mesoscopic shape transformation, similar to the one shown in Figure 2c along the direction 1-5. The OR plane is formed at ~3.5, which for typical LCs [6] corresponds to ~10 V/m.

Discussion
We studied numerically the annihilation of master surface enforced highly charged pair {defect,

Discussion
We studied numerically the annihilation of master surface enforced highly charged pair {defect, antidefect}, bearing topological charges {3/2, −3/2}. A Landau-de Gennes mesoscopic approach in terms of the tensor nematic order parameter was used. We assumed quasi-two-dimensional order, which in practice is sensible for relatively thin nematic cells. We chose geometrical setups and boundary conditions that could be realized experimentally in future work. For example, highly charged defect structures could be enforced using, e.g., the Atomic Force Measurement scribing method [30], or plasmonic photoalignment technique [32]. In simulations we enforce a pair of |m| = 3/2 defects via the command surface scribed pattern. We used anchoring strength coefficients that are strong enough to prevent annihilation of the pair, however weak enough to allow decomposition [33] of each surface pattern favored |m| = 3/2 defect into three daughter elementary defects bearing |m| = 1/2. Furthermore, we confined the pair in a rectangular confinement region. One would expect mutual annihilation of the two daughter assemblies on varying relevant control parameters. However, the chosen setup enabled a nontrivial annihilation scenario of TDs, also involving creation of additional defects, absorption of TDs within edges of the confining substrate, and one dimensional "explosion" of defect cores. We illustrate two qualitatively different annihilation channels, either by progressively weakening the relative importance of the surface anchoring field or increasing external in-plane electric field. Both processes could be realized experimentally by using relatively simple experimental setups. In our setting, the reference structures exhibit eight elementary TDs, where two additional defects are introduced by the rectangular confinement. In the anchoring weakening-driven annihilation, only four daughter defects were mutually annihilated, and the remaining two were annihilated with confinement enabled TDs. On the contrary, in E-driven annihilation only one pair of daughter defects were mutually annihilated and the remaining four defects were annihilated via TDs, which were created via creation and depinning of additional pairs {defect, antidefect}.
Note that in realistic LC samples the characteristic nematic elastic constants are different from each other. Nematic elasticity is commonly expressed in terms of splay (K 11 ), twist (K 22 ), bend (K 33 ), and saddle-splay (K 24 ) elastic constants, introduced in the Frank-Oseen uniaxial description [34]. In our geometry the elastic distortions involving twist and saddle-splay deformations are absent. Therefore, only K 11 and K 33 are expected to make an impact on defect structures. In conventional LCs it holds K 33 > K 11 . In this case the defect rearrangement would be carried out via LC structural variations that prefer the splay nematic deformation. Furthermore, this elastic anisotropy would break the mirror symmetry along the y-axis observed in Figure 4a. Namely, the core structures of m = 1/2 and |m| = −1/2 defects would be different, because the relative contributions of bend and splay deformations are different in these defects.
On approaching a second order N-SmA phase transition K 33 diverges, whereas K 11 remains finite. Because any arrangement of TDs in our study involves bend elastic distortions, TDs would annihilate on increasing K 33 . This elastic anisotropy would, for sure, at least quantitatively influence the annihilation scenario. The defect rearrangement would be realized via LC rearrangements preferring the splay deformation. However, for our scenario we do not expect a qualitatively different rearrangement of TDs from those presented in Figures 4 and 6.
Nevertheless, it is a demanding task to obtain a numerically "realistic" result of elastic anisotropy on rearrangement of TDs. In this paper we work in the approximation of equal elastic constants, represented by a temperature independent constant L. In reality, several different elastic constants would appear, which determine the relative importance of different symmetry-allowed elastic contributions expressed in terms of Q and its spatial derivatives. For example, in Cartesian coordinates (x,y,z) = (x 1 , x 2 , x 3 ), the expansion up to the 2nd order in Q is commonly expressed as [35,36] f e = f (2) e = L 1 Q ij,k Q ij,k + L 2 Q ij,j Q ik,k + L 3 Q ij,k Q ik,j + L 4 Q ij Q ij,kk + L 5 Q jk Q ik,jk .
Here we use the summation convention over repeated indices, commas indicate spatial derivatives (e.g., Q ij,k = ∂Q ij ∂x k ), and the superscript (2) in f (2) e indicates the 2 nd order expansion contributions in the elastic free energy density. Far from defect cores, the constants L 1 -L 5 could be mapped to the more familiar Frank elastic constants. However, for the expansion up to the 2 nd order in Q it holds that K 11 = K 33 [35,36], while in most LCs these constants can be significantly different. Higher order terms need to be added [35,36] to lift this degeneracy, for example f e = f (2) e + L 6 Q kl Q ij,l Q ij,k . However, according to recent studies in lyotropic chromonic LCs [24], even taking into account the above mentioned elastic anisotropy, one could not exactly reproduce experimentally the observed core structures of |m| = 1/2 TDs. Therefore, current theoretical modeling does not correctly describe cores of single |m| = 1/2 TDs, and the core structure details might influence annihilation scenarios. For this reason, we performed our study at the proof-of-principle level.
In our selected geometry pairs of defects {−1/2,1/2} are typically created as an intermediate step and they actively mediate external stimulus driven annihilation processes. To estimate the free energy costs of a pair {−1/2,1/2} creation we use the Lyuksyutov constraint [26,37] in which the cubic term in f c is neglected. This is sensible approximation for LCs exhibiting weakly first order I-N phase transition. With this in mind, it follows f c ∼ a 0 (T − T * )trQ 2 + c trQ 2 2 , where f c is minimized in the nematic phase for trQ 2 = a 0 (T * −T) 2c , corresponding to f c ∼ − a 2 0 (T−T * ) 2 4c . The main penalty of introducing a defect is "melting" of LC order within the core of TDs. Note that for a pair {−1/2,1/2}, the nematic director in the far field is spatially homogeneous because the total charge of the pair equals to zero. The free energy cost of introducing a single |m| = 1/2 defect of length h is roughly πξ 2 b h, where ∆V ∼ πξ 2 b h estimates the defect core volume. For a typical nematic LC it holds that a 0 ∼ 0.05 10 6 J Km 3 , c ∼ 10 6 J m 3 , ξ b ∼ 30 nm [23]. Therefore, deep in the nematic phase it costs roughly 2∆F/(k b T * ) ∼ 10 3 to create a pair of defects, where we take |T − T * | = 10 K, h ∼ ξ b , T * ∼ 330 K, and k b is the Boltzmann constant.

Conclusions
In summary, we studied the external stimulus-enforced annihilation of highly-charged assembly of command surface enforced TDs. In bulk, these TDs would simply mutually annihilate on varying an appropriate control parameter. However, we showed that in appropriate confinement these processes could be extremely complex and display qualitatively different realizations. In particular, we showed that when the characteristic confinement size and separation of highly charged defects are comparable, additional defect pairs are created and involved in the annihilation processes. In these events the total topological charge of the system is conserved in accordance with topological charge conservation laws. Our simulations exhibit a rich diversity of defect configurations for which only a relatively small set of control parameters was varied.