Properties of twisted topological defects in 2D nematic liquid crystals

Abstract

Topological defects are one of the most conspicuous features of liquid crystals. In two dimensional nematics, they have been shown to behave effectively as particles with both charge and orientation, which dictate their interactions. Here, we study “twisted” defects that have a radially dependent orientation. We find that twist can be partially relaxed through the creation and annihilation of defect pairs. By solving the equations for defect motion and calculating the forces on defects, we identify four distinct elements that govern the relative relaxational motion of interacting topological defects, namely attraction, repulsion, co-rotation and co-translation. The interaction of these effects can lead to intricate defect trajectories, which can be controlled by setting relevant timescales.

---

I. Introduction

Topological defects are of great interest in many areas of physics, with prominent examples in cosmology, astrophysics, and crystal structures.¹ Their existence can be hugely influential on the behaviour of a system, notably for phase transitions in condensed matter.² Beyond the physical context, topological defects are now also recognized to be an important aspect of morphogenesis in biological systems.^3–7 Topological defects are a prominent feature of liquid crystals, fluid materials in which the components have a broken rotational symmetry.^8–10 Indeed their identification led to the discovery of nematic phases, where the broken rotational symmetry is captured by a director field. Furthermore, topological defects can be used to determine the material properties of nematic materials.^11–15

In two-dimensional nematic liquid crystals, topological defects are point-like disclinations^1,9,10 with a well-defined half-integer charge and orientation.^16,17 These defects have an energetic cost as they disturb the order of the liquid crystal resulting in a Coulomb-like interaction potential between defects.^1,9,10 This results in opposite signed topological defects being attracted to each other and eventually annihilating. The situation is more interesting in “active” liquid crystals, where the continuous insertion of mechanical stress keeps the system out of equilibrium.^18,19 Active liquid crystals can exhibit a state referred to as “active turbulence” featuring the spontaneous generation of flows along with the proliferation of topological defects.²⁰ These active turbulent states generally achieve a steady state density of defects when the elastic and active forces balance.^21–24 Furthermore, the interaction between defects has been reported to lead to defect ordered states.^25–33

Elastic interactions between defects are present in both active and passive nematics, and have been shown to play an important role in the torques defects exert on each other.^16,17,33 Vromans & Giomi¹⁶ introduced an efficient method for identifying the orientation of half-integer topological defects in nematic liquid crystals, along with describing their interaction energy, results which have since been expanded upon by Tang & Selinger.¹⁷ In particular, these studies identify an elastic energy dependence on the relative orientation of topological defects, which also results in defects taking curved paths when repelling or attracting one another.^16,17

In this manuscript we investigate the behaviour of twisted topological defects in two-dimensional nematics. We call a topological defect twisted if the corresponding director field depends on the distance from the defect centre. After characterising individual twisted defects, we extend our analysis to the interaction between two twisted defects. Importantly, the corresponding elastic energy is not necessarily periodic as a function of the relative defect orientation. Instead, it can increase as the square of the relative defect orientation as defects are rotated relative to each other. Furthermore, we identify four distinct behaviours in the motion of pairs of twisted defects, namely the combination of attraction and repulsion with co-translation and co-rotation. We show that the origin of these behaviours lies in the combination of the signs of the defect charges and their twist relative to the background nematic texture. These features can lead to complex relaxation dynamics that can be controlled by temperature as we illustrate with an example of four defects.

II. Defects out of equilibrium

A liquid crystal can be described by the director field n = (cos(θ), sin(θ)). The director denotes the local average orientation of the anisotropic molecules, where θ is the angle relative to some background frame of reference. In the case of a nematic, the molecules are elongated and identical under reversal, thus any physics must be invariant under the transformation n → −n.

The energy of a nematic liquid crystal depends on variations of the director in space, which in two dimensions can take the form of bend or splay, each with their own energetic cost. For simplicity, we will consider the limit where the corresponding elastic constants are equal. Then, the total energy of a nematic texture is given by the Frank free energy^1,9,10

where K is the single elastic constant.

Topological defects in two dimensional nematics are points, such that the director winds by a multiple of π on a closed path encircling them, i.e..^1,9,10 Here is referred to as the topological charge of the defect. At the center (or core) of a defect, the director is not well defined. For an isolated defect with charge k, the director field minimising the Frank free energy is

θ = kϕ + θ₀. (2)

Here ϕ is the polar angle around the core and the phase, θ₀, is a constant that depends on the reference frame. Note that θ is independent of the distance r from the defect core. In the case of half integer charge, k = ±1/2, the defect is not rotationally symmetric and the phase θ₀ is related to the defect orientation.^16,17

For a configuration containing multiple defects, the director field minimising the Frank free energy is given by

where ϕ_i is the polar angle around the core of defect i with charge k_i and Θ₀ is the global phase.¹ The single global phase controls the phase at the core of every defect simultaneously, thus eqn (3) only describes the nematic texture around a set of defects which are in phase with each other. While this is usually sufficient in describing passive liquid crystals, in which the nematic quickly adopts the minimum energy configuration where defects are in phase with each other, this is not necessarily the case in out of equilibrium scenarios, such as those commonly found in active nematics.

In active nematics, there is a constant injection of mechanical stress at the microscopic scale which leads to a proliferation of defects, see Fig. 1a. Using the positions of the defects observed in Fig. 1a we construct a nematic texture using eqn (3). Here we see that the single fitting parameter, Θ₀, is not sufficient to match the observed orientation of every defect, see Fig. 1b. This is because the active forces drive defects from their equilibrium positions and orientations, leading to many defects out of phase with each other.

Fig. 1 Defects in a nematic out of thermodynamic equilibrium. (a) Snapshot of a microtubule-based active nematic featuring multiple positive (thick red dashes) and negative (blue triads) topological defects. From ref. 33. (b) Equilibrium nematic texture generated around the positions of the defects in the active nematic by eqn (2). (c) Nematic texture generated by eqn (5), which recreates the relative phases of each of the defects.

In order to capture the situation depicted in Fig. 1a, we introduce the idea of a twisted topological defect, in which the phase is a function of the radius r. Its general form is given by

θ = kϕ + τf(r) + θ₀, (4)

where f(r) is a monotonic function. Without loss of generality, we set ∂_rf(r) < 0, and f(0) = 1. The twist amplitude τ controls the defect orientation at its core relative to the nematic texture at r → ∞ rather than to any external frame of reference. It should be noted here that the introduced twist is not the same as the conventional liquid crystalline twist, n·(∇ × n), found in cholesteric phases, which by construction is zero in two dimensions.

We can describe the nematic texture around a set of twisted defects by writing

where ∂_rf_i(r) < 0, , and f_i(0) = 1 for all i. With this description, we can control the phase of each defect individually by adjusting the values of τ_i and f_i(r). To analyze the nematic order field of Fig. 1a, we use a single functional form of f_i for all i to obtain a unique set of values {τ_i}. Explicitly, f_i(r) = exp(−10r/L), where L is the size of the quadratic domain. Introducing both the positions and the orientations of the defects in Fig. 1a, into eqn (5) generates the nematic texture given in Fig. 1c, which clearly differs from the one in Fig. 1b and is much closer to the experimental data.

III. Single twisted defect

In order to explore the elastic effects of the defect twists introduced in eqn (4), we shall look at simple configurations of defects with high values of twist, starting with a single defect.

The nematic texture around a single twisted defect does not minimise the Frank free energy, leading to an increase in the energy by compared to the defect (2). Therefore, if the nematic texture around a defect is allowed to relax, in the long time limit f(r) = const. However, when fixing the phase along two boundaries at different radii, see Fig. 2a (blue lines), the twist might not fully relax, resulting in a frustrated configuration, see Movie S1 (ESI†).

Fig. 2 Single twisted defects. (a) Director field (thin red lines) and corresponding Schlieren texture (grey) of a +1/2 defect given by eqn (4) and (8) with r₀/L = 5/512, r₁ = 250/512, and τ = π. Thick red dashes indicate the defect orientation at the fixed boundaries (blue dashed lines). (b) Landau-De Gennes energy E_LdG of nematic textures initialised as in (a) as a function of the twist amplitude τ and the defect core radius ε. (c) E_LdG as a function of τ for two values of ε. Dotted lines indicate E ∼ τ². (d) Schematic of the unzipping removing a π twist between the boundaries. The defects nucleate close to the inner fixed boundary (blue dashed lines) and follow pink lines removing π twist and annihilating.

When exploring the behaviour of nematics numerically, it is advantageous to work with the nematic tensor Q_ij = S(n_in_j − δ_ij/2), where S is the order parameter, instead of the director n.^9,10 In contrast to the director n, the tensor Q is well defined even at the defect core, where the order parameter S = 0. Then, the Frank energy (1) is replaced by the Landau–De Gennes energy

The Landau–De Gennes energy introduces an additional parameter ε, often referred to as the defect core radius. The parameter ε describes the length scale over with the order parameter S varies.

The relaxation dynamics of the nematic can then be described by

where we have introduced the rotational viscosity γ. We numerically integrate eqn (7) using a finite difference scheme on an N × N grid with N = 512. When simulating stationary defects, we employ fixed boundary conditions, , indicated by blue dashed lines in Fig. 2(a and d) and 3(a, b and e). At the edges of the simulation we impose Neumann boundary conditions, which leave the director free to rotate at the boundary. In this case, there are no topological restrictions on the net defect charge in the simulation domain as is imposed by the Poincaré–Hopf theorem in the case of fixed boundary conditions.^1,16 In all cases studied in this work, the simulation domains are sufficiently large such that defects have at least a distance L/4 to the boundaries. In the case of repelling defects, which would eventually approach a boundary, simulations are stopped, when the distance of a defect to the nearest boundary is smaller than L/4. We scale all lengths by the size L and time by K/γL². All energies are presented in units of the elastic constant K. In keeping with previous works on the topic,^16,17 we have opted to focus on the elastic interactions between defects and to neglect hydrodynamic effects introduced by backflow or active stresses.

Fig. 3 Pairs of pinned twisted defects. (a) Director field (thin red lines) and corresponding Schlieren texture (grey) of a +1/2 and a −1/2 defect with twists τ₁ = −τ₂ = π, see text for details. Thick red dash and blue tripod indicate the defect orientations at the fixed boundaries (blue dashed lines). The defect cores are at ±1/8 along the x-axis, ε= 0.001, and r₀/L = 5/512. (b) Same as (a), but for two +1/2 defects. (c) Landau–De Gennes energy E_LdG of nematic textures initialised as in (a and b) as a function of the twist amplitude Δτ = τ₁ − τ₂ and defect core radius ε. (d) E_LdG as a function of Δτ for two values of ε. Dotted lines indicate E ∼ Δτ². (e) Schematic of the unzipping process corresponding to a ±1/2 defect pair. A pair of defects is nucleated close to each fixed boundary (blue dashed lines) and follow the pink lines, removing 2π of twist before annihilating.

We start by examining the relaxation of a twisted nematic liquid crystal between two fixed boundaries at r = r₀ and r = r₁. The nematic texture was initialised by imposing a director field according to eqn (4) with a linear twist

f(r) = τ(r₁ − r)/(r₁ − r₀), (8)

giving the configuration shown in Fig. 2a. This texture is then allowed to relax according to eqn (7) to a stable configuration. Although many of the nematic configurations studied here may appear extremely twisted, it is worth remembering that all results are scaled by the size of the domain. This is because only the net twist between two boundaries is important. Thus the tightness of the twist and therefore the energy density, can be decreased by simply increasing the distance between the boundaries. Integer defects with a high degree of twist have been generated in experimental systems by applying a rotating magnetic field to the nematic liquid crystal 5CB³⁴ or by exposing nematic colloidal suspensions to a chemical gradient.³⁵

In Fig. 2b, we present the value of E_LdG after relaxation as a function of the defect core radius ε and the initially applied twist amplitude τ. As τ is increased, the energy exhibits distinct discontinuities at locations which depend on ε. This is highlighted in Fig. 2c displaying the Landau–de Gennes energy of the twisted nematic as a function of τ for two values of ε. The energy increases as ΔE ∼ τ² before reaching a critical value, at which point it drops.

These drops are associated with the creation and annihilation of a defect pair that is able to “unzip” some of the twist between the two boundaries. Since the orientation at the boundaries is fixed, the difference in orientation between a point at r = r₀ and a point at r = r₁ has a fixed value. In contrast, in the course of the relaxation process, the twist between these points can change by increments of π. This is a consequence of the periodic nature of the nematic director.

For example, consider the nematic texture shown in Fig. 2a. Along a straight radial path from r = r₀ to r₁, the difference in twist is Δθ₀ = τ[f(r₁) − f(r₀)] = −τ. If the twist is increased such that τ → τ + π, the orientations at the boundary have not changed, whereas the nematic texture now stores an additional twist of π along the path. This extra π twist can only be released by nucleating a pair of oppositely charged defects, which during the relaxation process trace a path around the central boundary to annihilate with each other, see Fig. 2d and Movie S2 (ESI†). This effect has been observed previously for integer defects of liquid crystals in a quasi 2D “dowser” state.³⁴ There is an energy barrier for nucleating a defect pair. The height of this barrier depends on the defect core radius ε. Only if the energy stored in the twisted nematic is above this nucleation energy, can twist be released through this mechanism.

The critical value of τ tells us the maximum twist a nematic texture can elastically tolerate for a given scale. For example, the configuration shown in Fig. 2a has a twist of τ = π between two boundaries separated by a distance of r₁ − r₀ = 245/512. From Fig. 2b we see this is stable for values of ε ≲ 0.003. In the microtubule based active nematic, the defect core radius is estimated to be ε ∼ 10 μm,³³ hence a configuration similar to that shown in Fig. 2a would only be stable for systems larger than L ∼ 1.6 mm. For the nematic liquid crystal 5CB, integer defects with high degrees of twist have been observed at length scales L ∼ 10 mm.³⁴

We conclude that when the defect core radius is very small compared to the distance between the two boundaries, it is possible for the nematic texture to store a stable twist greater than π.

IV. Defect interaction energy

We now turn our attention to the interaction of defects. Let us consider two fixed circular boundaries with radius r₀, each surrounding the core of a defect, and let the angle θ be given by

θ = k₁ϕ₁ + k₂ϕ₂. (9)

Here, k_i, i = 1, 2, are the charges of the two defects and ϕ_i are the polar coordinates relative to the respective defects. This creates a nematic texture, which minimises the Frank free energy. Note that the value of θ is constant along the straight line connecting the two defect cores. We then add a linear twist to each defect according to eqn (8) with r₁ being half the inter-defect distance and f(r) = 0 for r > r₁. This configuration then relaxes according to eqn (7), while the director is fixed on the boundaries. Examples for the resulting nematic textures are shown in Fig. 3a, b and Movie S3 (ESI†).

In the following, we consider defects with k = ±1/2 and distinguish between defect pairs of opposite, Fig. 3a, or the same charge, Fig. 3b. In both cases, the defects have equal and opposite twists with respect to the background nematic, τ₁ = −τ₂. The apparent difference in chirality between the Schlieren textures of these two configurations comes from the fact that τ > 0 will rotate the Schlieren texture of opposite charges in opposite directions.

Two observations should be noted. First, the defects must have an applied twist of opposite sign, τ₁ = −τ₂. If the defects are twisted with the same sign they remain in phase with each other. In this situation, any twist introduced can be negated by a global rotation of the background nematic, therefore this type of twist is not preserved by the boundaries, see Movie S4 (ESI†). Secondly, in the long time limit only the relative applied twist of the defects is important. Defects can exchange twist by a global rotation of the background nematic and the lowest energy configuration is always that in which the twist is spread evenly between the two boundaries, see Movie S5 (ESI†). Hence, a configuration in which one defect is twisted by 2π will relax to the configuration in which both defects are twisted by ±π.

As in the single defect case, the net difference between the angle θ for points on the two boundaries is fixed. This implies that the relative twist between the two defects can again exceed π. After relaxation, the energy is qualitatively similar to the single defect case, Fig. 3c. It again features drops in energy associated with an “unzipping” of twist by the introduction of pairs of defects, which move around the boundaries and annihilate. Note that the form of the upper line of energy discontinuities results from trapping a defect close to the fixed boundary. Also, the relationship ΔE ∼ Δτ² is preserved and holds up to the points at which unzipping becomes energetically favourable, Fig. 3d. The quadratic dependence on the difference in twist amplitudes is consistent with previous observations made by Vromans & Giomi.¹⁶

Again the critical value of τ gives us an estimate of the length scale at which a certain phase difference between a pair of defects is stable. The configurations in Fig. 3a, b show defects with Δτ = 2π, which are stable for values of ε ≲ 0.003. Thus to observe a similar configuration in a microtubule based active nematic, the defect separation would need to be of the order of 1 mm.

During unzipping, defects are nucleated close to the boundaries around the defect cores, where the elastic energy density is highest, Fig. 3e. As the twist and the elastic energy associated with this process is equally shared between the two defects, two pairs of defects are nucleated simultaneously during the relaxation process, one close to each boundary. Each new defect pair then encircles the nearby boundary and annihilates, see Movie S6 (ESI†). Compared to the single twisted defect, this process requires double the defect nucleation energy and unzips 2π of twist. If the separation of the two boundaries is much greater than the defect core radius ε, there is not sufficient energy stored in the nematic texture to nucleate defects and a difference in angles between the two boundaries can increase above 2π. Importantly, this implies that the interaction energy between two defects is not necessarily a periodic function of their relative orientation as has been previously predicted.¹⁶

V. Defect motion

When the fixed boundaries are removed and the nematic texture is allowed to fully relax to a homogenous state, a ±1/2 defect pair attracts and annihilates, whereas defects with the same sign repel each other. As has been previously reported, when the defects have different orientations, they typically do not travel along a straight line toward or away from each other during the relaxation process.^16,17 Instead, defects of opposite signs co-rotate while attracting, Fig. 4a, b and Movie S7 (ESI†). This is a direct consequence of the fact that in order to annihilate, two defects must be in phase with each other, i.e. the intermediate texture can be described by eqn (3). Thus when two out of phase defects approach each other, they must first remove any relative twist between them before annihilating.²² Conversely, defects of the same sign co-translate while repelling each other, Fig. 4c, d and Movie S8 (ESI†). Formally, these co-rotations and co-translations are expressed by a rotation of the vector connecting the defect cores and by a translation of their mean position, respectively.

Fig. 4 Dynamics of out of phase defect pairs for different core radii. (a) Defect trajectories for the annihilation of a ±1/2 defect pair initialised with opposite orientations for different values of the core radii ε. Inset: The decay of the twist amplitude difference Δτ is independent of ε. (b) Defect separation as a function of time for the defect pairs shown in (a). Inset: Average velocity of the defects as a function of the defect core radius. (c) Same as (a) but for a repelling pair of +1/2 defects. (d) Same as (b) but for the defect pairs in (c).

This behaviour does not change qualitatively with the defect core radius ε. Also, the rate at which the twist amplitude difference Δτ is dissipated depends only weakly on the defect core radius ε, Fig. 4a and c (insets). In contrast, the translational speed of defects exhibits a more pronounced dependence on ε, Fig. 4b and d, and the velocity of the defects increases with the defect core radius, Fig. 4b and d (insets).

The weak or strong dependence on the defect core radius ε can be explained by the existence of two competing timescales that are relevant for the relaxation dynamics determined by eqn (7). The timescale associated with relaxation of the nematic director is given by T_θ = Δr²γ/K, where Δr is the typical length over which θ varies; in the present case the inter defect spacing. The timescale associated with relaxation of the order parameter is given by T_S = ε²γ/K. For the orientation difference between the defects to relax, only the director needs to vary, whereas the order parameter can remain stationary. Hence this process is dictated by T_θ and therefore independent of ε. In contrast, for a defect core to move, both the nematic director and the order parameter must change. This is because the order parameter drops from S ≈ 1 to S ≈ 0 as you approach the defect core. This process thus depends on T_S and hence ε.

As mentioned above, for fixed defects, a relative twist can be preserved only if the twists have opposite signs and we only discussed this case. Now that we consider the relaxation dynamics, however, also situations with equal signs of the twist should be analysed. This defines four distinct prototypic scenarios for interacting defect pairs: pairs with like or opposite charges and like or opposite twist amplitudes. The corresponding relaxation dynamics are displayed in Fig. 5 and Movies S7–S10 (ESI†). As above, the dynamics can be either rotational, Fig. 5a and d, or translational, Fig. 5b and c. This pattern is also apparent from the Schlieren textures of the initial nematic fields, which are either chiral for co-rotating configurations, Fig. 5a and d (insets), or achiral for co-translating configurations, Fig. 5b and c (insets). By comparing Fig. 4 and 5 it is clear that the degree of perpendicular motion of the defects depends more strongly on the difference in twist amplitudes Δτ than the size of the defect core radii ε.

Fig. 5 Dynamics of out of phase defect pairs for different initial twist amplitude differences following eqn (7). (a) Defects with opposite charges ±1/2 and equal twist amplitude τ; (b) defects with opposite charges ±1/2 and opposite twist amplitudes ±τ; (c) defects with equal charges 1/2 and equal twist amplitudes τ; (d) defects with equal charges and opposite twist amplitudes ±τ. Insets: Schlieren texture of the corresponding initial configurations. Parameter values: ε = 0.008 and initial τ = 0 (blue), 1 (violet), 2 (red), 3 (dark red), 4 (dark green), 5 (bright green). For repelling defects, trajectories are shown up to the point at which they approach the boundary closer than L/4 and the twist is almost entirely dissipated. From this point defects move directly away from each other.

VI. Forces on the defect cores

To better understand the motion of defects described in the previous section, we analyze in the following the forces and torques acting on twisted defect cores. In this case, the director field around a set of N twisted defects is conveniently written aswhere ϕ_i and r_i are the polar angle and distance relative to defect i, respectively, with i = 1,…,N. Note that this ansatz neglects a possible dependence of f_i on the polar angle ϕ_i. As above f_i is a monotonic function describing the local twist applied to defect i and τ_i are the twist amplitudes.

The interaction between such defects is governed by the Frank energy. It can be obtained by calculating the gradient of θ around each defect in the polar coordinates centred on that defect with basis set [[small phi, Greek, circumflex]_i, r̂_i], such that

We now focus on the case of two interacting defects. The square of the gradient of the angle θ can then be expressed as

Using the fact that the basis vectors are orthonormal and that [small phi, Greek, circumflex]₁[small phi, Greek, circumflex]₂ = r̂₁r̂₂ and [small phi, Greek, circumflex]₁r̂₂ = −r̂₁[small phi, Greek, circumflex]₂, this can be simplified to

The first four terms are the defect self-energies associated with the defects' topological charges and twists. The defect self-energy associated with twist is proportional to τ², as was observed in Fig. 2c and 3d. This was first observed by Vromans & Giomi¹⁶ and calculated exactly for a texture minimizing the Frank free energy by Tang & Selinger.¹⁷ The next terms describe the way the two defect charges and twists interact. The term proportional to [small phi, Greek, circumflex]₁[small phi, Greek, circumflex]₂ has a part proportional to the product of the defect charges that is akin to the Coulomb interaction between electrical charges. It is familiar from interacting defects in the absence of twist, τ = 0.^1,9,10 The contribution describing the interaction of the two twists is similar, however, the dependence on the radial coordinates is given by f′ instead of 1/r. The final term results from a coupling between twist and charge.

Let us examine the symmetries of the various interaction terms. To this end we assume without loss of generality that both defects are positioned at ±r along the x-axis. We will also assume that they have equal or opposite charges, k₁= ±k₂, as well as equal amplitudes, |τ₁| = |τ₂|, and distributions, f₁ = f₂. The term proportional to [small phi, Greek, circumflex]₁[small phi, Greek, circumflex]₂ is symmetric for x → −x and for y → −y. This then describes the energy density between the two defects and can lead to attraction and repulsion.

The final term describes the twist-charge interaction and is proportional to [small phi, Greek, circumflex]₁r̂₂, which is anti-symmetric in y, thus describes the relative energy density above and below the defects. In addition to this, the symmetries of the final term depend on the relative signs of the defect charges and twists. If the defects have either the same charge and same twist or opposite charge and opposite twist, it can be written as . This expression is anti-symmetric in both x and y and thus can lead to co-rotation of the defects. If in contrast the defects have either the same charge and opposite twist or the opposite charge and the same twist, the final term can be written as , which is symmetric in x and anti-symmetric in y, thus can lead to co-translation. This can be summarized by the rule defects will co-rotate if sgn(k₁k₂) = sgn(τ₁τ₂) and will co-translate if sgn(k₁k₂) = −sgn(τ₁τ₂). Furthermore, the direction of the co-translation or co-rotation depends on the sign of τ₁τ₂. These different symmetries induce the qualitatively distinct defect dynamics observed in the previous section.

To obtain the forces and torques on the defect cores, we differentiate the energy with respect to the defect positions x̲_i and twist amplitudes τ_i. We additionally assume here that the f_i remain fixed such that the time-dependence of all twist is captured by the coefficients τ_i. This allows us to write a set of over-damped dynamical equations for the positions of the defects and their twist amplitudes:

Here, μ and ν play the role of a translational and rotational mobility, respectively. We again scale all lengths by the size of the integration domain, L, and times by Kμ/L². The integrals are performed numerically on a 512 × 512 grid.

Let us consider the case , where r′ indicates the length scale over which the twist decays. We numerically integrate these equations to obtain the trajectories of defect pairs for different values of τ_i and k_i, Fig. 6. A comparison with Fig. 5 shows that our model captures the qualitative features of the dynamics of twisted defect pairs. It should be noted here that the dynamics cannot be expected to agree quantitatively as the form f_i(r_i) = exp(−r_i/r₀) is not a minimiser of the Frank free energy.

Fig. 6 Defect trajectories obtained from eqn (14) and (15) corresponding to Fig. 5. Parameter values: μ/ν = 1.95, initial τ = 0 (blue), 1 (violet), 2 (red), 3 (dark green), 4 (bright green).

VII. Multi defect dynamics

We now turn our attention to decay dynamics of more complex arrangements of twisted defects and their non-trivial trajectories. Consider an arrangement featuring a horizontal line of four equally spaced defects with charges −1/2, +1/2, +1/2, −1/2 as viewed from left to right. In a twist free texture, the defects would obviously annihilate as a pair on the left and a pair on the right. If we now apply twists with amplitudes +π/2, −π/2, −π/2, +π/2 this will result in the nematic texture shown in Fig. 7a. In this configuration, any two adjacent defects will co-rotate and adjacent pairs co-rotate in the opposite direction. When this texture relaxes, the defects take curved paths and annihilate as a pair on the left and a pair on the right, dashed lines Fig. 7c, d and Movie S11 (ESI†). If the twist of the central pair is now increased, such that the twists are now given by +π/2, −3π/2, −3π/2, +π/2, we obtain the texture given in Fig. 7b. In this configuration, the central pair of defects co-rotate to such a degree that they will switch places before annihilating with the initially furthest negative defect, solid lines Fig. 7c, d and Movie S12 (ESI†).

Fig. 7 Relaxation of two defect pairs. (a) Orientation field (red dashes) and corresponding Schlieren texture (grey) given by eqn (10) for four defects with charges −1/2, +1/2, +1/2, and −1/2 from left to right. The corresponding twist amplitudes are +π/2, −π/2, −π/2 and +π/2 with and r′ = 0.1. (b) Same as (a) but with twists +π/2, −3π/2, −3π/2 and +π/2. (c) Defect trajectories corresponding to nematic texture shown in (a) dashed lines and (b) solid lines from eqn (14) and (15). Parameter values μ/ν = 1.3. (d) Defect trajectories corresponding to (a) dashed lines and (b) solid lines from eqn (7). Parameter values: ε = 0.016. (e) Defect trajectories corresponding to nematic texture shown in (b) from eqn (14) and (15). Parameter values: μ/ν = 0.65 (dashed lines) and 0.97 (solid lines). (f) Defect trajectories corresponding to nematic texture shown in (b) from eqn (7). Parameter values: ε = 0.004 (dashed lines) and ε = 0.008 (solid lines).

Since the elastic energy of a twisted defect E ∼ τ², the additional twist of the central pair of defects will initially decay faster than that of the outer pair of defects. Depending on the values of the mobilities μ and ν, the defects will either reorient quickly or slowly compared to their translational motion. Hence by changing the relative mobilities, we can control the degree of co-rotation that occurs. When the rotational dynamics are fast, the higher twists decay quickly and the defects annihilate as a pair on the left and a pair on the right, dashed lines Fig. 7e. However, when the rotational dynamics are slowed down, the additional twist on the central +1/2 defects again causes them to switch places and annihilate with the initially furthest −1/2 defects, solid lines Fig. 7e.

The behaviour of the effective defect dynamics is again paralleled by the solutions to the relaxation dynamics given by eqn (7), Fig. 7f and Movies S12 and S13 (ESI†). There, the translational mobility is analogous to decreasing the defect core radius ε. Depending on its value, either adjacent pairs annihilate (dashed lines) or the central +1/2 defects move away from the nearest −1/2 defects and annihilate with the originally more distant −1/2 defects (solid lines). Excitingly, this means that the trajectories of relaxing defects can be controlled by adjusting the defect core radius, which in an experimental setting is linked to the temperature of the liquid crystal.

VIII. Conclusions

Topological defects in two dimensional nematics are often described as particles that exist within the liquid crystal and interact through their charge. More recently, this idea has been expanded to include the fact that defects have both charge and orientation.^16,17 The defect charges interact through a Coulomb like interaction, with like charges repelling and opposite charges attracting. Relative defect orientations have been shown to increase the interaction energy of topological defects and cause transverse motion as defects attract or repel.^16,17,33

In this paper, we have considered the situation where defects are not in phase with each other. In that case, which is common in active nematics far from equilibrium, the nematic orientations corresponding to the defects are incompatible with each other and the intermediate texture does not minimise the Frank free energy. In this case, we found that the nematic order can be captured by a superposition of defects with a radially dependent phase, which takes the form of a local twist to the nematic texture. We have also shown that the local orientation of a defect with respect to an external reference frame is typically not sufficient to characterise the relaxation dynamics of defects. Consequently, the relaxation dynamics depend on global features of the orientation field. For a single defect with fixed orientation angle at the boundaries, twist may persist, because the configuration has insufficient energy to generate a defect pair, which is necessary to “unzip” a π rotation of the director. These findings extend to situations with several defects.

In cases, where the orientation angle is not constrained at the boundaries and the nematic texture can relax to a defect-free state, the relaxation dynamics depend on global properties of the orientation field as well as on several time scales. The relative twist and charge interact and cause characteristic off-axis motion. The time scales are associated with changes in orientation that keep the order parameter constant and with defect core displacements, which necessarily involve changes in the order parameter. They are in turn set by two length scales, namely, the characteristic distance between defects and the defect core radius.

The qualitative dependence of defect motion on the core radius could play a significant role in the dynamics of active nematics. Indeed, active nematics are often associated with a constant density of defects and states with long range defect order have been speculated upon both experimentally and theoretically.^25–33 The process of spontaneous defect ordering is likely to be strongly affected by the two time scales identified above.

Conflicts of interest

There are no conflicts to declare.

Appendix

A Fitting defect twists to experimental data

In the following, we describe how the twist of topological defects can be determined from experimental data. The positions and orientations of defects are identified using the methodology outlined by Vromans et al.¹⁶ These are then converted into a position, r̲_i, charge, k_i, and phase, ψ_i for each defect i. In the neighbourhood of defect i, the nematic texture can then be approximated by

θ = k_iϕ_i + ψ_i. (A1)

If we equate the nematic texture in the neighbourhood of defect j with the general, twisted nematic texture given in eqn (5) we can write

We now introduce the matrices, ϕ and F, that will allow us to compute the twists τ_i. The components of ϕ are ϕ_ij = arctan2(ŷ·(r̲_j − r̲_i), x̂·(r̲_j − r̲_i)) and yield the polar angle of defect j around defect i, for i ≠ j. Furthermore, we set ϕ_ii = 0. In turn, the components of F are F_ij = f_i(r_ij), where r_ij is the distance between the cores of defect i and defect j, with F_ii = 1 as per the boundary conditions introduced in the main text.

Finally, we define δ_i as the deviation of the local phase at defect i from the equilibrium

Since F is symmetric, this equation can easily be inverted provided a suitable choice of f_i(r). Thus the twists required to reproduce a set of defects with given positions, charges and phases are given by

Since we introduced the same number of constraints, ψ_i, as parameters, τ_i, we are left with a single parameter, Θ₀, which can be used to fit the remaining background nematic.

B Supplementary movie captions

• Supp_Mov_1.mp4 (ESI†): relaxation of a single +1/2 defect between fixed boundaries with a twist given by eqn (8) and τ = 2π. This is the same as is shown in Fig. 2a.

• Supp_Mov_2.mp4 (ESI†): animated schematic of unzipping of twist around a +1/2 defect.

• Supp_Mov_3.mp4 (ESI†): relaxation of twist between two boundaries each containing an oppositely charged defect. Each defect is initialised with a linear twist of π with opposite sign, resulting in a net twist between the boundaries of Δτ = 2π. This twist is topologically preserved.

• Supp_Mov_4.mp4 (ESI†): relaxation of twist between two boundaries each containing an oppositely charged defect. Each defect is initialised with a linear twist of π with same sign, resulting in a net twist between the boundaries of Δτ = 0. This twist is not topologically preserved, and is removed by a global phase change.

• Supp_Mov_5.mp4 (ESI†): relaxation of twist between two boundaries each containing an oppositely charged defect. The right hand defect is initialised with a linear twist of π, resulting in a net twist between the boundaries of Δτ = π. This twist is topologically preserved, but shares equally between the boundaries, leading to each defect having a twist of τ = ±π/2.

• Supp_Mov_6.mp4 (ESI†): animated schematic of unzipping of twist around a defect pair.

• Supp_Mov_7.mp4 (ESI†): relaxation and annihilation of a pair of oppositely charged twisted defects. The defects are twisted out of phase, i.e. τ₁ = −τ₂. The pair co-rotate before annihilating.

• Supp_Mov_8.mp4 (ESI†): relaxation and repulsion of a pair of twisted +1/2 defects. The defects are twisted out of phase, i.e. τ₁ = −τ₂. The pair co-translate while repelling.

• Supp_Mov_9.mp4 (ESI†): relaxation and annihilation of a pair of oppositely charged twisted defects. The defects are twisted while remaining phase, i.e. τ₁ = τ₂. The pair co-translate before annihilating.

• Supp_Mov_10.mp4 (ESI†): relaxation and repulsion of a pair of twisted +1/2 defects. The defects are twisted while remaining in phase, i.e. τ₁ = τ₂. The pair co-rotate while repelling.

• Supp_Mov_11.mp4 (ESI†): annihilation of two pairs of oppositely charged defects. From left to right, the charges are −1/2, +1/2, +1/2, −1/2, and the initial twists are τ = +π/2, τ = −π/2, τ = −π/2, τ = +π/2. Each defect annihilates with the closest oppositely charged defect. This is the same data as shown in Fig. 7a and d.

• Supp_Mov_12.mp4 (ESI†): annihilation of two pairs of oppositely charged defects. From left to right, the charges are −1/2, +1/2, +1/2, −1/2, and the initial twists are τ = +π/2, τ = −3π/2, τ = −3π/2, τ = +π/2. The defect core radius is small, ε = 0.004, each defect annihilates with the closest oppositely charged defect. This is the same data as shown in Fig. 7b and f.

• Supp_Mov_13.mp4 (ESI†): annihilation of two pairs of oppositely charged defects. From left to right, the charges are −1/2, +1/2, +1/2, −1/2, and the initial twists are τ = +π/2, τ = −3π/2, τ = −3π/2, τ = +π/2. The defect core radius is large, ε = 0.008, each defect annihilates with the closest oppositely charged defect. This is the same data as shown in Fig. 7b and f.

Acknowledgements

This work was funded by the SNSF through the NCCR for Chemical Biology. We thank Jyothishraj Nambisan for providing us the experimental image of the active nematic shown in Fig. 1a and Pau Guillamat for insightful discussions.

References

1. P. Chaikin and T. C. Lubensky, Principles of condensed matter physics, Cambridge University Press, Cambridge, 1995.
2. J. M. Kosterlitz and D. J. Thouless, J. Phys. C: Solid State Phys., 1973, 6, 1181.
3. A. Singer, L. Boucheron, S. H. Dietze, K. E. Jensen, D. Vine, I. McNulty, E. R. Dufresne, R. O. Prum, S. G. J. Mochrie and O. G. Shpyrko, Sci. Adv., 2016, 2, E1600149.
4. T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Nature, 2017, 544, 212.
5. K. Kawaguchi, R. Kageyama and M. Sano, Nature, 2017, 545, 327.
6. Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun and K. Keren, Nat. Phys., 2021, 17, 251.
7. P. Guillamat, C. Blanch-Mercader, K. Kruse and A. Roux, 2020, bioRxiv 2020.06.02.129262.
8. F. C. Frank, Faraday Discuss., 1958, 25, 19.
9. P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, Clarendon Press, Oxford, 1995.
10. M. Kleman and O. D. Lavrentovich, Soft Matter Physics: An Introduction, Springer, New York, 2003.
11. S. D. Hudson and E. L. Thomas, Phys. Rev. Lett., 1989, 62, 1993.
12. J. Brugues, J. Ignés-Mullol, J. Casademunt and F. Sagués, Phys. Rev. Lett., 2008, 100, 037801.
13. R. Zhang, N. Kumar, J. L. Ross, M. L. Gardel and J. J. De Pablo, Proc. Natl. Acad. Sci. U. S. A., 2017, 115, E124.
14. C. Blanch-Mercader, P. Guillamat, A. Roux and K. Kruse, Phys. Rev. Lett., 2021, 126, 028101.
15. C. Blanch-Mercader, P. Guillamat, A. Roux and K. Kruse, Phys. Rev. E, 2021, 103, 012405.
16. A. Vromans and L. Giomi, Soft Matter, 2016, 12, 6490.
17. X. Tang and J. V. Selinger, Soft Matter, 2017, 13, 5481.
18. S. Ramaswamy, Annu. Rev. Condens. Matter Phys., 2010, 1, 323.
19. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143.
20. T. Sanchez, D. N. Chen, S. J. DeCamp, M. Heymann and Z. Dogic, Nature, 2012, 491, 431.
21. L. Giomi, M. J. Bowick, X. Ma and M. C. Marchetti, Phys. Rev. Lett., 2013, 110, 228101.
22. L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. C. Marchetti, Philos. Trans. R. Soc., A, 2014, 372, 20130365.
23. L. Giomi, Phys. Rev. X, 2015, 5, 031003.
24. P. Guillamat, J. Ignés-Mullol and F. Sagués, Nat. Commun., 2017, 8, 564.
25. S. J. DeCamp, G. S. Redner, A. Baskaran, M. F. Hagan and Z. Dogic, Nat. Mater., 2015, 14, 1110.
26. E. Putzig, G. S. Redner, A. Baskaran and A. Baskaran, Soft Matter, 2016, 12, 3854.
27. A. Doostmohammadi, M. F. Adamer, S. P. Thampi and J. M. Yeomans, Nat. Commun., 2016, 7, 10557.
28. A. U. Oza and J. Dunkel, New J. Phys., 2016, 18, 093006.
29. D. J. G. Pearce, Phys. Rev. Lett., 2019, 122, 227801.
30. S. Shankar and M. C. Marchetti, Phys. Rev. X, 2019, 9, 041047.
31. K. Thijssen, L. Metselaar, J. M. Yeomans and A. Doostmohammadi, Soft Matter, 2020, 16, 2065.
32. D. J. G. Pearce, New J. Phys., 2020, 22, 063051.
33. D. J. G. Pearce, J. Nambisan, P. W. Ellis, A. Fernandez-Nieves and L. Giomi, 2020, arXiv:2004.13704.
34. P. Pieranski, S. Čopar, M. H. Godinho and M. Dazza, Eur. Phys. J. E: Soft Matter Biol. Phys., 2016, 39, 121.
35. S. Hernàndez-Navarro, P. Tierno, J. A. Farrera, J. Ignés-Mullol and F. Sagués, Angew. Chem., Int. Ed., 2014, 53, 10696.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1sm00825k

---