Orientation of topological defects in 2D nematic liquid crystals

TL;DR

This paper introduces a tensor-based formalism to describe the orientation of topological defects in 2D nematic liquid crystals, providing new insights into defect interactions and dynamics.

Contribution

It proposes an alternative tensor formalism for defect orientation, improving upon previous vector-based models, and applies it to analyze defect interactions and motion.

Findings

Tensor formalism effectively describes defect orientations.

Orientation influences defect interaction energies.

Numerical simulations show defect motion depends on orientation.

Abstract

Topological defects are an essential part of the structure and dynamics of all liquid crystals, and they are particularly important in experiments and simulations on active liquid crystals. In a recent paper, Vromans and Giomi [Soft Matter, 2016, 12, 6490] pointed out that topological defects are not point-like objects but actually have orientational properties, which strongly affect the energetics and motion of the defects. That paper developed a mathematical formalism which describes the orientational properties as vectors. Here, we agree with the basic concept of defect orientation, but we suggest an alternative mathematical formalism. We represent the defect orientation by a tensor, with a rank that depends on the topological charge: rank 1 for a charge of +1/2, rank 3 for a charge of -1/2. Using this tensor formalism, we calculate the orientation-dependent interaction between defects, and we present numerical simulations of defect motion.

Equations

$$F=\frac{K}{2}\int d^{2}r|\nabla \theta {|}^2,$$

$$\theta =k\varphi +{\theta }_0,$$

$$\theta =\varphi (\text{mod }\pi ).$$

$$\psi =\frac{{\theta }_0}{1-k}\left(\text{mod }\frac{\pi }{|1-k|}\right).$$

$$\mathbf{p}=\frac{\nabla \cdot (\mathbf{nn})}{|\nabla \cdot (\mathbf{nn})|}.$$

$$n_i{n}_j=\frac{{\delta }_{ij}}{2}+\frac{r_i{p}_j+r_j{p}_i-(\mathbf{r}\cdot \mathbf{p}){\delta }_{ij}}{2|\mathbf{r}|}.$$

$$\nabla \cdot (\mathbf{nn})=\mathbf{n}(\nabla \cdot \mathbf{n})+(\mathbf{n}\cdot \nabla )\mathbf{n}=\mathbf{n}(\nabla \cdot \mathbf{n})-\mathbf{n}\times (\nabla \times \mathbf{n}).$$

$$T_{ijk}=\frac{2}{3}\left(p_i^{(1)}{p}_j^{(1)}{p}_k^{(1)}+p_i^{(2)}{p}_j^{(2)}{p}_k^{(2)}+p_i^{(3)}{p}_j^{(3)}{p}_k^{(3)}\right).$$

$$|\mathbf{T}{|}^2\equiv T_{ijk}{T}_{ijk}=1.$$

$$n_i{n}_j=\frac{{\delta }_{ij}}{2}+\frac{T_{ijk}{r}_k}{|\mathbf{r}|}.$$

$${\partial }_k(n_i{n}_j)=\frac{1}{|\mathbf{r}|}\left[T_{ijk}-\frac{T_{ijl}{r}_l{r}_k}{|\mathbf{r}{|}^2}\right].$$

$$⟨{\partial }_k(n_i{n}_j)⟩=\frac{1}{2|\mathbf{r}|}{T}_{ijk}.$$

$$\mathbf{T}=\frac{⟨\nabla (\mathbf{nn})⟩}{|⟨\nabla (\mathbf{nn})⟩|}.$$

$$T_{ijk}=$$

$$\frac{⟨{\partial }_i(n_j{n}_k)+{\partial }_j(n_k{n}_i)+{\partial }_k(n_i{n}_j)⟩}{\sqrt{⟨{\partial }_a(n_b{n}_c)+{\partial }_b(n_c{n}_a)+{\partial }_c(n_a{n}_b)⟩⟨{\partial }_a(n_b{n}_c)+{\partial }_b(n_c{n}_a)+{\partial }_c(n_a{n}_b)⟩}},$$

$$n_x{n}_x$$

$$=\frac{1+\cos (2k\varphi )\cos (2(1-k)\psi )-\sin (2k\varphi )\sin (2(1-k)\psi )}{2},$$

$$n_x{n}_y$$

$$=\frac{\sin (2k\varphi )\cos (2(1-k)\psi )+\cos (2k\varphi )\sin (2(1-k)\psi )}{2}.$$

$$\theta =\varphi (\text{mod }2\pi /n).$$

$$\psi =\frac{{\theta }_0}{1-k}\left(\text{mod }\frac{2\pi }{n|1-k|}\right).$$

$$T_{ij}=\frac{1}{\sqrt{2}}\left(p_i^{(1)}{p}_j^{(1)}+p_i^{(2)}{p}_j^{(2)}-{\delta }_{ij}\right),$$

$$\mathbf{n}=\frac{\sqrt{2}\mathbf{T}\cdot \mathbf{r}}{|\mathbf{r}|},$$

$${\partial }_j{n}_i=\frac{\sqrt{2}}{|\mathbf{r}|}\left[T_{ij}-\frac{T_{ik}{r}_k{r}_j}{|\mathbf{r}{|}^2}\right].$$

$$⟨{\partial }_j{n}_i⟩=\frac{1}{\sqrt{2}|\mathbf{r}|}{T}_{ij}.$$

$$\mathbf{T}=\frac{⟨\nabla \mathbf{n}⟩}{|⟨\nabla \mathbf{n}⟩|},$$

$$T_{ij}=\frac{⟨{\partial }_i{n}_j+{\partial }_j{n}_i⟩}{\sqrt{⟨{\partial }_a{n}_b+{\partial }_b{n}_a⟩⟨{\partial }_a{n}_b+{\partial }_b{n}_a⟩}}.$$

$$\theta (\mathbf{r})\approx \{\begin{array}{ll}-\frac{1}{2}{\tan }^{-1}(y/x)+{\theta }_1& \text{for }\mathbf{r}\text{ near }(0,0),\\ +\frac{1}{2}{\tan }^{-1}(y/(x-R))+{\theta }_2& \text{for }\mathbf{r}\text{ near }(R,0).\end{array}$$

$${\nabla }^2\theta =0.$$

$$w=\frac{z-\gamma r_{\text{core}}}{\gamma z-r_{\text{core}}},z=\frac{r_{\text{core}}(w-\gamma )}{\gamma w-1},$$

$$\gamma =\frac{R+\sqrt{R^2-4r_{\text{core}}^2}}{2r_{\text{core}}}.$$

$$w_{\text{min}}=\frac{R^2-2r_{\text{core}}^2-R\sqrt{R^2-4r_{\text{core}}^2}}{2r_{\text{core}}^2},$$

Full text

\makeFNbottom

[TABLE]

††footnotetext: a* Liquid Crystal Institute, Kent State University, Kent, OH 44242, USA; E-mail: jselinge@kent.edu*

1 Introduction

Topological defects are common in many areas of physics, from high-energy physics and cosmology to crystal structure and superconductivity.1, 2 Indeed, the importance of topological defects was recognized by the 2016 Nobel Prize in Physics. In the science of liquid crystals, topological defects are often used to recognize phases, and they determine the structure and dynamics of many liquid-crystal phases with technological applications.3

In recent years, topological defects have particularly been studied in the context of two-dimensional (2D) active nematic liquid crystals,4 which typically occur in systems of rod-like macromolecules on a surface. This research has identified two new features of topological defects. First, in active systems, defects of topological charge $+1/2$ are constantly moving, with motion driven by the activity of the material. By contrast, defects of topological charge $-1/2$ are almost at rest, with only slow diffusive motion. Second, both $+1/2$ and $-1/2$ defects have characteristic orientations. For $+1/2$ defects, the orientation can been seen in the comet-like texture of the director field, and it is also the direction of the driven motion. For $-1/2$ defects, the orientation can be seen in the three-fold symmetric texture of the director field around the defect. Surprisingly, experiments and simulations have both shown that active systems can have long-range order in the orientation of the defects.5, 6, 7 The defects tend to maintain alignment along some spontaneously chosen axis, even as they are constantly created and annihilated.

In order to understand this long-range order, one must consider the general concept of defect orientation. A geometric question is: How can the orientation of a defect be described mathematically? Further physical questions are: How does the orientation affect the interaction between two defects, or the motion of a defect? Although these questions are motivated by studies of active nematic systems, they are relevant to all nematic liquid crystals, even equilibrium phases.

In the context of 3D nematic liquid crystals, Čopar et al.8, 9 investigate defect orientation by defining the splay-bend parameter. This parameter provides an excellent method to visualize the orientation of a disclination line, but it does not give a mathematical expression for the orientation. More recently, in the context of 2D nematic liquid crystals, Vromans and Giomi10 argue that the orientational properties of nematic defects can be described by vectors, and give an explicit expression for the vectors. Based on the vector construction, they calculate the orientation-dependent interaction between defects. Furthermore, they model the relaxational dynamics of interacting defects, and find that the trajectory depends strongly on the defect orientation.

The purpose of our current work is to examine the concept of defect orientation in 2D nematic liquid crystals in more detail. Through this study, we partially agree and partially disagree with the work of Vromans and Giomi. First, in Sec. 2, we show that their vector formalism is quite reasonable for $+1/2$ defects. However, for other defect charges, we represent the defect orientation by a tensor, with a rank that depends on the topological charge of the defect and on the symmetry of the underlying phase. In a 2D nematic phase, a $+1/2$ defect is represented by a tensor of rank 1 (i.e. a vector), while a $-1/2$ defect is represented by a tensor of rank 3. For a general $n$-atic phase (with an orientational order parameter of $n$-fold symmetry), a topological defect of charge $k$ is represented by a tensor of rank $n|1-k|$.

Next, in Sec. 3, we investigate the interaction between neighboring defects. Using conformal mapping, we determine the director field around two defects with arbitrary orientations, and calculate the elastic energy associated with that director field. The result is different from the interaction reported by Vromans and Giomi. In Sec. 4, we construct partial differential equations to model the relaxational dynamics of a 2D nematic phase, and use these equations to simulate the annihilation of $\pm 1/2$ defects with arbitrary initial orientations. For the dynamics, we agree with the results of Vromans and Giomi: The defect trajectories are quite similar to those reported in their paper, and these trajectories depend sensitively on defect orientation. Finally, in Sec. 5, we discuss defect orientation as a concept for understanding the physics of 2D nematic liquid crystals.

2 Orientation and tensor structure of defects

Like Vromans and Giomi, we consider a 2D nematic liquid crystal. This phase has orientational order of the molecules along the local director $\mathbf{n}(\mathbf{r})=(\cos\theta(\mathbf{r}),\sin\theta(\mathbf{r}))$. Because the orientational order is two-fold symmetric, the director $\mathbf{n}$ is equivalent to $-\mathbf{n}$, and hence $\theta$ is only defined modulo $\pi$. For that reason, the topological defects in this phase are disclination points, about which the director rotates through a multiple of $\pi$, so that $\oint d\theta=2\pi k$. Here, $k$ is the topological charge of the defect, which must be an integer or half-integer, positive or negative.

The elastic free energy of the 2D nematic phase is the Frank free energy. In the approximation of equal Frank elastic constants, this free energy can be written as

[TABLE]

where $K$ is the single Frank constant. The local minimum of $F$ corresponding to a defect of topological charge $k$ at the origin is given by

[TABLE]

where $\phi=\tan^{-1}(y/x)$ is the angle in polar coordinates. The angle $\theta_{0}$ represents an arbitrary overall rotation of the director about the $z$-axis. Figure 1 shows several examples of these defects, with $\theta_{0}=\pi/3$ in all cases.

To characterize the orientation of a defect, we can ask: Where does the director $\mathbf{n}$ point radially outward from (or inward toward) the defect? This occurs when the angle $\theta$ satisfies

[TABLE]

Solving Eqs. (2) and (3) simultaneously, we find that these special radial directions are given by $\theta=\phi=\psi$, where

[TABLE]

Hence, the radial directions can be described by the vector $\mathbf{p}=(\cos\psi,\sin\psi)$, which is defined up to rotations through $\pi/|1-k|$. This is precisely the defect orientation vector defined by Vromans and Giomi. In Fig. 1, the red arrows show the rotationally equivalent $\mathbf{p}$ vectors for each topological charge $k$.

We now must consider how the $\mathbf{p}$ vector is related to the director field around the defect, for different values of $k$.

2.1 Topological charge $k=+1/2$

For topological charge $k=+1/2$, Eq. (4) shows that the angle $\psi$ is defined modulo $2\pi$. As a result, the defect orientation vector $\mathbf{p}=(\cos\psi,\sin\psi)$ is a single-valued vector, as shown in Fig. 1(d). It should be possible to determine $\mathbf{p}$ from the director field around the defect, and conversely, to determine the director field around the defect from $\mathbf{p}$.

To determine $\mathbf{p}$ from the director field, Vromans and Giomi show that

[TABLE]

To determine the director field around the defect, we construct a covariant expression for the tensor $\mathbf{nn}$ in terms of the defect orientation vector $\mathbf{p}$ and the position $\mathbf{r}$,

[TABLE]

An explicit calculation shows that this expression gives the same director field as Eq. (2) with $k=1/2$. As a specific check, we can see the behavior in the $\mathbf{p}$ direction out from the defect: If $\mathbf{r}/|\mathbf{r}|=\mathbf{p}$, then $n_{i}n_{j}=p_{i}p_{j}$, so that the director field points radially outward (or inward). As another check, the divergence of Eq. (6) is $\partial_{i}(n_{i}n_{j})=p_{j}/(2|\mathbf{r}|)$, and hence the vector $\mathbf{p}$ is the normalized version of $\nabla\cdot(\mathbf{nn})$, consistent with Eq. (5). Hence, this covariant expression provides a way to work with the director field in terms of $\mathbf{p}$ and $\mathbf{r}$.

Because the covariant expression of Eq. (6) is equivalent to Eq. (2) with $k=1/2$, it is valid in the same regime: outside a core radius $r_{\text{core}}$, up to the distance where the director field is affected by different mechanisms beyond the central defect, such as other defects or boundaries. These equations are derived with the approximation of equal Frank elastic constants. If the Frank constants are unequal, then the director field would have corrections with higher powers of $\mathbf{r}/|\mathbf{r}|$.

As Vromans and Giomi point out, one physical interpretation of Eq. (5) is in the context of active nematic liquid crystals: The stress tensor includes an active term proportional to $\mathbf{nn}$, with positive or negative sign depending on whether the material is contractile or extensile. The active force on the defect is then in the direction $\nabla\cdot(\mathbf{nn})$. Hence, $\mathbf{p}$ gives the direction of defect motion in an active nematic.

We can suggest another physical interpretation in the context of flexoelectricity.11, 12 The divergence term can be expanded as

[TABLE]

The right side of this equation is the standard expression for the flexoelectric polarization of a liquid crystal, in the case where the splay and bend flexoelectric coefficients $e_{1}$ and $e_{3}$ are equal (which is an approximation similar to the approximation of equal Frank constants). Hence, a $+1/2$ defect is a point of highly concentrated flexoelectric polarization, in the direction given by $\mathbf{p}$. As discussed by Čopar et al.,8, 9 the splay-bend parameter is the divergence of this flexoelectric polarization.

2.2 Topological charge $k=-1/2$

For topological charge $k=-1/2$, Eq. (4) shows that the angle $\psi$ is defined modulo $2\pi/3$. As a result, the defect orientation vector $\mathbf{p}=(\cos\psi,\sin\psi)$ is really a triple-valued vector, as shown in Fig. 1(b). In other words, three vectors are equivalent representations of the same defect orientation: $\mathbf{p}^{(1)}$, $\mathbf{p}^{(2)}$, and $\mathbf{p}^{(3)}$, which are related to each other by rotations through $2\pi/3$. This triple-valued vector is not necessarily a problem; physicists often need to work with mathematical objects that are multiple-valued. (Indeed, the director field $\mathbf{n}$ is an example of a double-valued vector field, because $\mathbf{n}$ and $-\mathbf{n}$ are equivalent representations of the same orientational distribution.)

Vromans and Giomi extract $\mathbf{p}$ from the director field around a $-1/2$ defect through a two-step derivation: They first define an intermediate vector $\mathbf{p}^{\prime}$, which depends on the arbitrary choice of coordinate system, and then use $\mathbf{p}^{\prime}$ to calculate $\mathbf{p}$, which does not depend on the coordinate system. Here, we propose an alternative approach, which is consistent with their formalism but may be clearer because it explicitly uses the symmetry of the system.

Our alternative approach is based on a higher-rank tensor. In general, an object with $r$-fold symmetry can be represented by a single-valued, completely symmetric tensor of rank $r$. Hence, for a $-1/2$ defect, we consider the tensor $\mathbf{T}$ with components

[TABLE]

This tensor is invariant under rotations of $2\pi/3$, and it is normalized so that

[TABLE]

To determine the director field around the defect, we construct a covariant expression for the tensor $\mathbf{nn}$ in terms of the defect orientation tensor $\mathbf{T}$ and the position $\mathbf{r}$,

[TABLE]

An explicit calculation shows that this expression gives the same director field as Eq. (2) with $k=-1/2$. As a specific check, we can see the behavior in the $\mathbf{p}^{(1)}$ direction out from the defect: If $\mathbf{r}/|\mathbf{r}|=\mathbf{p}^{(1)}$, then $n_{i}n_{j}=p^{(1)}_{i}p^{(1)}_{j}$, so that the director field points radially outward (or inward). The same is true for $\mathbf{p}^{(2)}$ and $\mathbf{p}^{(3)}$.

To determine the defect orientation tensor from the director field, we take the gradient of Eq. (10),

[TABLE]

By averaging this expression in a region around the defect, so that $\langle r_{l}r_{k}/|\mathbf{r}|^{2}\rangle=\delta_{lk}/2$, we obtain

[TABLE]

Hence, the defect orientation tensor can be determined from the director field as the normalized average

[TABLE]

A more symmetric version of this expression is

[TABLE]

Equations (13) and (14) are equivalent for the exact director field of Eq. (10), but Eq. (14) might be more suitable for numerical calculations in general. Once we have $T_{ijk}$, we can determine the three directions $\mathbf{p}^{(1)}$, $\mathbf{p}^{(2)}$, and $\mathbf{p}^{(3)}$ through the following construction: Define a test unit vector $\mathbf{b}=(\cos\beta,\sin\beta)$, calculate the scalar $f(\beta)=T_{ijk}b_{i}b_{j}b_{k}=\frac{1}{2}\cos 3(\beta-\psi)$, and find the maxima of $f(\beta)$. Those maxima occur at the angles $\beta=\psi$ (mod $2\pi/3$). This construction identifies the principal orientation of the third-rank tensor $T_{ijk}$, just as an eigenvector calculation identifies the principal axes of a second-rank tensor. In that sense, it can be considered as a generalized eigenvector calculation.

We have tested this construction on sample textures containing $-1/2$ defects, and confirm that it identifies the three directions where the director field points outward (or inward). Hence, it identifies the orientation of the defects, in a way that does not depend on any choice of coordinate system.

2.3 Other topological charges in a 2D nematic phase

For a defect with an arbitrary topological charge $k$, Eq. (4) shows that the angle $\psi$ is defined up to rotations through $\pi/|1-k|$. Hence, the defect has $2|1-k|$-fold rotational symmetry, with an orientation represented equivalently by $2|1-k|$ distinct $\mathbf{p}$ vectors. Instead of this multiple-valued vector, we can also describe it by a single-valued, completely symmetric tensor of rank $2|1-k|$. This tensor can be constructed from the $\mathbf{p}$ vectors by generalizing Eq. (8).

The general relationship between the director field and the defect orientation tensor can be seen by expanding the tensor $\mathbf{nn}$. From Eqs. (2) and (4), an explicit calculation gives

[TABLE]

Here, each factor of $\cos(2k\phi)$ or $\sin(2k\phi)$ can be expressed in terms of $2|k|$ factors of $\mathbf{r}/|\mathbf{r}|$. Likewise, each factor of $\cos(2{(1-k)}\psi)$ or $\sin(2{(1-k)}\psi)$ can be expressed in terms of $2|1-k|$ factors of $\mathbf{p}$, or equivalently in terms of a defect orientation tensor with rank ${2|1-k|}$. A covariant, single-valued description of the defect orientation requires a tensor of that rank.

One important special case is a defect of topological charge $k=+1$, as shown in Fig. 1(e). In that case, the tensor rank is $2|1-k|=0$, and hence the tensor is just a scalar. This result is physically reasonable, because the defect is an isotropic object, which has no special directions going outward from the core. Like other defect charges, the $+1$ defect has a parameter $\theta_{0}$ in the director field. However, the significance of this parameter is different for a $+1$ defect than for other defect charges. For a $+1$ defect, the parameter $\theta_{0}$ determines whether the director field points radially (splay deformation), tangentially (bend deformation), or somewhere in between. This is a scalar property of the defect, which does not change as one moves around the defect core. It is not related to an orientation of the defect.

2.4 Generalization to $n$-atic phases

Apart from the 2D nematic phase, researchers often consider phases with other types of orientational order. If the orientational order parameter has $n$-fold symmetry, then the phase is called $n$-atic. The case $n=1$ is a polar phase, such as a ferromagnet or ferroelectric, with a vector order parameter. The case $n=2$ is a nematic phase, as discussed above. The case $n=6$ is a hexatic phase, which commonly arises from 2D bond-orientational order.13

We can generalize the theory presented in this section to an $n$-atic phase with arbitrary $n$. The elastic free energy still has the form of Eq. (1), and the orientational order around a defect is still given by Eq. (2). For general $n$, the topological charge $k$ must be a positive or negative integer multiple of $1/n$, so that the director rotates through an angle of $2\pi/n$ around a loop about the defect. The special directions where the orientational order matches the outward radial alignment are now given by

[TABLE]

Solving Eqs. (2) and (17) simultaneously, we find that these special radial directions occur at $\theta=\phi=\psi$, where

[TABLE]

Once again, we can construct the vector $\mathbf{p}=(\cos\psi,\sin\psi)$, which is now defined up to rotations through $2\pi/(n|1-k|)$. Hence, there are $n|1-k|$ rotationally equivalent $\mathbf{p}$ vectors; i.e. the defect is an object with $n|1-k|$-fold rotational symmetry. It can therefore be represented by a completely symmetric tensor of rank $n|1-k|$.

Figure 2 shows the example of a polar phase, which is $n$-atic with $n=1$. The arrows can represent the magnetization in a ferromagnet, the electrostatic polarization in a ferroelectric, or the tilt in a smectic-C liquid crystal. Here, the topological charge $k$ must be an integer. For $k=+1$, the defect is characterized by a scalar (just as it is for the nematic case, or for any arbitrary $n$). Indeed, the example of Fig. 2(c) shows that this defect appears the same in all directions.

For $k=-1$, the defect has 2-fold rotational symmetry, as indicated by the double-headed red arrow in Fig. 2(a). The defect orientation is characterized by a tensor $\mathbf{T}$ of rank 2, with components

[TABLE]

with the $\delta_{ij}$ subtracted in order to make the tensor traceless. This tensor is invariant under rotations of $\pi$, and it is normalized so that $|\mathbf{T}|^{2}\equiv T_{ij}T_{ij}=1$. To determine the vector field $\mathbf{n}$ around the defect, we construct the covariant expression

[TABLE]

which is equivalent to Eq. (2) with $k=-1$. To determine $\mathbf{T}$ from $\mathbf{n}$, we take the gradient of Eq. (20),

[TABLE]

By averaging this expression in a region around the defect, so that $\langle r_{k}r_{j}/|\mathbf{r}|^{2}\rangle=\delta_{kj}/2$, we obtain

[TABLE]

Hence, $\mathbf{T}$ is the normalized average

[TABLE]

or more symmetrically,

[TABLE]

The $\mathbf{p}$ vectors are $\pm 1$ times the eigenvector of $\mathbf{T}$ with positive eigenvalue.

3 Interaction of defects

In this section, we determine the director field around two defects of arbitrary orientation in a 2D nematic phase. We then use that director field to calculate the interaction energy as a function of the distance and relative orientation between the defects.

Suppose we have a defect of charge $-1/2$ at the origin, and a defect of charge $+1/2$ at the position $(R,0)$, assuming $R>0$. We would like to specify the orientations of these two defects independently, so that

[TABLE]

Unfortunately, we cannot just add these two solutions to obtain the solution for $\theta(\mathbf{r})$ everywhere. This sum would have defects with the correct charges in the correct locations, but it would not have the correct defect orientations. Rather, the constant terms $\theta_{1}$ and $\theta_{2}$ would just combine to give an overall constant, which would not allow us to fix the relative orientation of the two defects.

To find the director field around these defects, we must solve a differential equation with the appropriate boundary conditions. The Euler-Lagrange equation associated with the Frank free energy (1) is just Laplace’s equation

[TABLE]

Each defect has some core radius $r_{\text{core}}$, such that the nematic order is disrupted inside the core. We apply boundary conditions at the core radius $r_{\text{core}}$ around each of the defects. Hence, we must solve Laplace’s equation on the full $(x,y)$ plane except the two defect cores.

For this solution, we use the technique of conformal mapping, as illustrated in Fig. 3. We make a conformal transformation from the complex plane of $z=x+iy$ to the complex plane of $w=u+iv$, where

[TABLE]

with

[TABLE]

With this transformation, the circular boundary of radius $r_{\text{core}}$ about $z=0$ maps onto a circular boundary of radius $w_{\text{max}}=1$ about $w=0$, and the circular boundary of radius $r_{\text{core}}$ about $z=R$ maps onto a circular boundary of radius

[TABLE]

also about $w=0$. Hence, we must solve Laplace’s equation between two concentric circles in the complex $w$ plane. The solution is

[TABLE]

where

[TABLE]

Transforming back into the $(x,y)$ plane, this solution becomes

[TABLE]

In the limit of small core radius $r_{\text{core}}$, it can be approximated by

[TABLE]

Here, the inverse tangents and the additive constant are the usual expression for the director field around two defects. The term proportional to $\delta\theta$ is a new term, which is required to specify the relative orientation of the two defects.

Figure 4 shows the director field of Eq. (33), in the case where $r_{\text{core}}/R=0.1$. The $-1/2$ defect on the left has the fixed orientation $\theta_{1}=\pi/2$, which implies $\psi_{1}=\pi/3$ by the argument in the previous section. The $+1/2$ defect on the right has an orientation that rotates from $\theta_{2}=0$ to $\pi$, which implies $\psi_{2}=0$ to $2\pi$. In Fig. 4(a), the defects clearly have the optimal relative orientation, and the director field has the usual form as the sum of inverse tangents. As the $+1/2$ defect rotates, the director field becomes more distorted. When the defect rotates through a full circle, the texture does not return to its original form, because extra distortion has wound up throughout the director field.

Note that we can specify the defect orientations at the specific core radius $r_{\text{core}}$. We cannot specify the orientations at the centers of the defects, because the problem becomes mathematically undefined in the limit of $r_{\text{core}}\to 0$.

We can generalize the form of Eq. (33) to describe arbitrary defect charges $k_{1}$ and $k_{2}$ at arbitrary positions $\mathbf{R}_{1}=(x_{1},y_{1})$ and $\mathbf{R}_{1}=(x_{2},y_{2})$. This generalization gives

[TABLE]

where

[TABLE]

The visualization of the director field shows the same type of behavior as in Fig. 4. There is an optimal director field when the relative orientation has $\delta\theta=0$, and the configuration becomes more distorted as $\delta\theta$ increases.

To calculate the elastic free energy associated with the distorted director field, we put the expression for $\theta(\mathbf{r})$ from Eq. (34) into the Frank free energy of Eq. (1), and integrate over the plane out to the system size of $R_{\text{max}}$. The integral is done in Mathematica, using Cartesian coordinates such that $x_{1}=-\frac{1}{2}R$, $x_{2}=\frac{1}{2}R$, and $y_{1}=y_{2}=0$, over the domain $-R_{\text{max}}<x<-\frac{1}{2}R-a$, $-\frac{1}{2}R+a<x<\frac{1}{2}R-a$, and $\frac{1}{2}R+a<x<R_{\text{max}}$, with $-\infty<y<\infty$. The result is

[TABLE]

Here, the first term is the usual energy cost of a defect pair with net topological charge $(k_{1}+k_{2})$, which diverges logarithmically with system size $R_{\text{max}}$ unless the net charge is zero. The second term is the usual Coulomb-like logarithmic interaction between two defects, which is repulsive for like charges and attractive for opposite charges. The third term is a new contribution, which favors orientational alignment between the defects toward the optimal orientation of $\delta\theta=0$. It creates an aligning torque

[TABLE]

and this torque decreases as $\pi K/\log(|\mathbf{R}_{1}-\mathbf{R}_{2}|/r_{\text{core}})$ when the defect separation is much greater than the core radius.

In Eq. (36), the orientational interaction is expressed in terms of $\delta\theta$ rather than in terms of the defect orientation vectors $\mathbf{p}$ or tensors $\mathbf{T}$ defined in the previous section. This is necessary because the interaction is not a single-valued function of $\mathbf{p}$ or $\mathbf{T}$. As an example, the textures in Figs. 4(a) and 4(d) have the same $\mathbf{p}$ and $\mathbf{T}$, but clearly the texture in Fig. 4(d) is more distorted and has a higher elastic free energy. When $\delta\theta$ is small, it may be possible to express the orientational interaction in terms of $\mathbf{p}$ or $\mathbf{T}$. However, that cannot work when $\delta\theta$ is large and the texture is wound up, as in Fig. 4(d).

Our defect interaction of Eq. (36) can be compared with the work of Vromans and Giomi. 10 They calculate the defect interaction in two ways. In their first method, they construct a field $\theta(\mathbf{r})$ that linearly interpolates between the arctangents around the two defects, and then calculate the free energy associated with this field. For two $+1/2$ disclinations separated by a distance $d$ in the $x$ direction, in a square $L\times L$ domain, they find an orientation-dependent part of the free energy that scales as $(KL/d)(\delta\psi-\pi)^{2}$ (where $\delta\psi-\pi$ is $2\delta\theta$ in our notation). One should note that their linear interpolation is not a minimizer of the free energy, and hence the free energy that they calculate is higher than our free energy, with a different dependence on system size and defect separation. We would argue that the interaction between defects is only defined when the director field between the defects is a minimizer of the free energy (or perhaps has Casimir-like fluctuations about the minimizer). To our understanding, the free energy that they calculate is the free energy of a particular choice of director field, but it cannot be considered as a property of the defects.

In their second method, Vromans and Giomi use an image construction to model like-sign defects with charge $k_{1}=k_{2}\equiv k$, in the limit of large system size $L$ and large defect separation $d$, and find the orientation-dependent part of the interaction as $F=-\pi Kk^{2}\log(1-\mathbf{p}_{1}\cdot\mathbf{p}_{2})$. Although this result is expressed in terms of the vectors $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$, it really only applies over a limited domain of $0<\delta\psi<2\pi$ around the minimum; it does not describe arbitrary windings of $\mathbf{p}_{1}$ or $\mathbf{p}_{2}$ through a full circle. The quadratic minimum of this function is similar to our defect interaction, but it shows deviations from quadratic behavior that we do not find, and it does not show the logarithmic decay of our interaction.

We have done numerical simulations to check the result of Eq. (36). In these simulations, we construct a 2D hexagonal lattice of spins interacting through the energy $E=J\sum_{\langle i,j\rangle}{[1-\cos(2(\theta_{i}-\theta_{j}))]}$, which is a discretized approximation to the Frank free energy of Eq. (1). We fix the position and orientations of two defects by placing effective “particles” on certain plaquettes between lattice sites. Each particle has strong anchoring to fix the spins on its boundary in a defect configuration, with topological charge of $\pm 1/2$, and with specified orientation. We then use a relaxation method to minimize the total energy over the spins on all of the non-anchored sites. Through this method, we find the minimum total energy as a function of relative orientation. The results are consistent with the predicted quadratic dependence on $\delta\theta$. We were not able to check the distance dependence of the orientational interaction, because the logarithmic decay is very slow in comparison with accessible length scales.

4 Motion of defects

In this section, we investigate the motion of two opposite-charged defects as they annihilate each other, to determine how the motion depends on defect orientation.

In order to describe a system in which the defects are free to move, we generalize the theory to represent the magnitude and direction of nematic order by a tensor

[TABLE]

In that formalism, the free energy can be expressed as

[TABLE]

Away from defects, where gradients of $\theta$ are small, the bulk value of the scalar order parameter is $S=\sqrt{a/b}$. At the defect points, where $\theta$ is singular, the scalar order parameter goes to $S=0$. In each core around a defect point, the scalar order parameter varies over a length scale $r_{\text{core}}=\sqrt{L/a}$.

To model the time evolution of nematic order, we use the equations for pure relaxational dynamics

[TABLE]

where $\gamma_{1}$ is the rotational viscosity, and $\mathbf{Q}$ is constrained to be symmetric and traceless. For the initial condition, we use a director field containing two defects, with topological charges $k_{1}=-1/2$ and $k_{2}=+1/2$ and arbitrary initial orientations $\theta_{1}$ and $\theta_{2}$, as found in Eq. (34). We assume that the initial $S(\mathbf{r},t)$ has the bulk value everywhere except in the defect cores, with the functional form

[TABLE]

This expression is physically motivated, in that it goes to the bulk value of $S=\sqrt{a/b}$ away from the defects, and it goes to zero linearly at the defects, but the exact form is arbitrary. We use open boundary conditions, at which $\mathbf{Q}$ is free and the normal derivatives vanish. We solve the differential equations in Mathematica, iterating forward in time until the defects annihilate each other. At each time, we find the defect positions by searching for points where $S(\mathbf{r},t)$ vanishes, and then find the defect orientations by the procedure described in Sec. 2.

Figure 5 shows a series of snapshots of the time evolution of the $\mathbf{Q}$ tensor field. The system begins with the two defects at an unfavorable relative orientation, with $\delta\theta=\pi$. In the early stage of the dynamic process, the defects rapidly rotate into the optimal relative orientation. While they rotate, they also move in the $y$ direction, which is transverse to their separation in the $x$ direction. Once they reach the optimal relative orientation, the dynamics becomes much slower. In this stage of the process, the defects move straight toward each other. After they annihilate each other, they leave a defect-free configuration, which eventually becomes uniform.

To demonstrate the influence of defect orientation, Fig. 6 shows the defect trajectories for several values of the initial $\delta\theta$. In this figure, the symbols represent the positions at equally spaced times, and hence the spacing between the symbols indicates the defect velocity. For initial $\delta\theta=0$, the defects are already at the optimal relative orientation at the beginning of the calculation. In that case, they move straight toward each other. The motion is initially rapid as the $\mathbf{Q}$ tensor field relaxes from a somewhat arbitrary initial configuration, then it slows down once the system reaches an almost-stable configuration with two defects, then it accelerates as the defects grow closer and the attractive force between them increases. By contrast, for initial $\delta\theta\not=0$, the initial stage of motion involves both rotation and translation in the $y$ direction, transverse to the inter-defect separation, until the optimal relative orientation is reached. For larger initial $\delta\theta$, the amount of translation is greater. Once the defects have the optimal relative orientation, their motion becomes much slower, and then later accelerates after the defects grow closer.

We emphasize that the initial motion is in the $y$ direction, transverse to the inter-defect separation, in spite of the fact that the defect interaction of Eq. (36) depends only on the magnitude of the separation. In other words, the defects move in the transverse direction although the interaction provides no force in the transverse direction. This behavior is an important feature of the motion of objects with internal orientation, which have an anisotropic drag as they move through a medium. In this sense, the motion of a defect is analogous to the motion of a sailboat, which can move transverse to the wind because of the orientation of the boat and the sail.

Our results for defect motion are actually quite similar to the results of Vromans and Giomi. 10 They also find curved trajectories, which are induced by the initial relative orientation of the defects. Their calculations of the dynamic evolution of the $\mathbf{Q}$ tensor are not affected by the issues involving the defect interaction discussed in the previous section.

5 Discussion

In this paper, we have examined the concept of defect orientation, which was initially developed by Vromans and Giomi. 10 Through this study, we partially agree and partially disagree with their work. We agree with them about the vector description of defects with topological charge $+1/2$, and about the motion of defects. We suggest that a tensor formalism provides a clearer way to describe defects with other topological charges, although it is consistent with their vector formalism. We disagree with them about the interaction between defects. Of course, despite these specific differences, we recognize their contribution of introducing this concept into the theoretical physics literature.

We must emphasize that defect orientation is not a topological invariant like defect charge. Indeed, it is not a topological concept at all. Rather, it is a geometric feature of defects, which will certainly change as a function of time. In that respect, it is analogous to defect position, which also changes as a function of time. Physicists often speak of defects as if they were effective “particles,” which can move around inside a liquid crystal. We argue that they should be considered as particles with orientation as well as position. They can exhibit both rotational and translational motion, and these two types of motion are coupled together.

We expect that the concept of defect orientation can be generalized in several ways. One important generalization will be to connect it back to active nematic liquid crystals. While this concept was inspired by experiments and simulations on active materials, our calculations have so far only considered the case of equilibrium liquid crystals. Once it is combined with theories of active nematics, there will certainly be a coupling between vector orientation of $+1/2$ defects and active motion, and there may also be new ways to understand the long-range ordering of defect orientation. A further generalization will be to 3D nematic liquid crystals. In general, 3D nematics have a more complex set of defects than 2D nematics, with both hedgehog points and disclination lines. These defects have their own types of orientation, as investigated by the graphical visualizations and topological arguments of Čopar et al.,8, 9 and we expect that these orientations can be understood through new tensor constructions. Finally, the concept of defect orientation can be connected with other aspects of liquid crystal theory, including backflow, interaction with colloidal particles, and background alignment of the director field, leading to new insights into defect behavior.

We would like to thank A. Baskaran and L. Giomi for helpful discussions. This work was supported by National Science Foundation Grant No. DMR-1409658.