At present, liquid crystals are widely used as displays (low-molar-mass liquid crystals) and high-strength engineering plastics (liquid crystalline polymers) by applying their anisotropy of material constants. On the other hand, because a liquid crystal has both solid and liquid properties, it has a possibly much wider range of applications, as do solid and liquid. Therefore, studies on new applications of liquid crystals are needed.
When a nematic liquid crystal is observed under a polarizing microscope, a region in which molecules orient parallel or perpendicular to the direction of polarization appears as a dark field, a region in which molecules orient ±45˚ to the direction of polarization appears as a bright field, and a region in which molecules orient in the other direction appears as a gray field. Figure 1(a) is a photograph of a nematic state right after the phase transition from an isotropic state of N-(p-methoxybenzylidene)-p’-butyl-aniline (MBBA), a typical low-molar-mass nematic liquid crystal. We observe the Schlieren texture peculiar to liquid crystals. There are some points in Figure 1(a) where black lines intersect with each other. These points are called defects   , around which spatial distortions of molecular orientation fields are generated. Figure 1(b) shows the representative defect orientation states. At the defect cores, the direction of orientation is discontinuous, and the orientation state is singular.
The defects often generated in manufacturing liquid crystal products cause degradation of the productivity and performance of the products. However, it is experimentally known that a pair of defects such as those in Figure 1(b), which have different molecular orientations, attract and finally annihilate each other  . This defect annihilation changes the orientation direction, so that a flow is expected to be generated  -  . Applying thermal energy induces defect generation because the orientation order at defects is lower. Therefore, generating defects artificially at arbitrary locations by applying thermal energy would inevitably generate flows there, leading to the development of a new type of microactuator that converts thermal energy into kinetic energy.
A theory that is able to describe the molecular orientation state satisfactorily should include the effects of long-range and short-range molecular order as well as the effect of flow. Two theories, the Leslie-Ericksen theory   and the Doi theory   , are well-known and have been widely used for predicting the dynamics of liquid crystals. From the aforementioned viewpoint, however, the Leslie-Ericksen theory, which does not include the short-range order effect, and the Doi theory, which does not account for the long-range order effect, are out as candidates. Marrucci and Greco  have expanded the Maier-Saupe potential  used in the Doi theory to include both the short-range and long-range
Figure 1. Liquid crystal defect structures. (a) A photograph of a nematic state of MBBA. (b) Molecular orientation configurations around defect cores: Ellipses represent the molecules, and small circles are the defect cores.
order effects and applied it to the original Doi theory  . Feng et al.  have proposed an expression for the stress tensor and to complete the theory.
Our objectives are to study the flow induced by the annihilation of a pair of defects and to estimate its magnitude using the aforementioned theory as a constitutive equation. A number of studies have simulated the liquid crystal flow during the annihilation of paired defects     . Experiments for the moving speed of defects under electric fields have also been performed   . However, they do not aim to develop actuators but focus on the movement of defects rather than the liquid crystal flow.
To calculate the orientational distribution functions (ODFs) at many points in a physical space, we have to account for both computational accuracy and time. In this paper, we: 1) approximate the ODF with a series of spherical harmonic functions, 2) study the minimum number of expanded terms required to simulate the orientation state properly, and 3) finally estimate the magnitude of the induced flow during the annihilation of a pair of defects to obtain useful data that can contribute to developing new actuators.
2.1. Governing Equations
The evolution equation for the ODF f is written as
Here, t is the time, k the Boltzmann constant, T the absolute temperature, u the unit vector parallel to a liquid crystal molecule, the differential operator on a unit sphere, and κ the velocity gradient tensor. and V are the rotary diffusivity and Marrucci-Greco nematic potential  expressed as
where D is the rotary diffusivity in an isotropic state, the differential operator in physical space, and S the order parameter tensor defined as
δ is the unit tensor.
The conservation equations for the isothermal slow flow of liquid crystals are the continuity and linear momentum equations
where v is the velocity vector, ρ is the fluid density, and p is the pressure. τ is the extra stress tensor derived by Feng et al.  expressed as
Here, c is the number density of liquid crystal molecules, U the dimensionless nematic potential intensity, li a parameter indicating the long-range order effect of molecules, ζ the drag coefficient of rotary molecules. A and Q are the second and fourth moments of the ODF f, respectively, expressed as
2.2. Computational Procedure and Boundary Conditions
Let us consider a two-dimensional square area with a side length of H, shown in Figure 2, and put a pair of defects with different orientational states on two points whose coordinates are P (0.3H, 0.5H) and Q (0.7H, 0.5H). For the coordinate system in Figure 2, the components of the velocity gradient tensor κ and order parameter tensor S are
where u and w are the x and z components of the velocity vector v.
For computation of the orientation field, we substitute Equation (10) and Equation (11) into Equation (1) to obtain the ODF f. In this study, we approximate f with a finite series of spherical harmonic functions Ylm(u)    .
Here, Clm(t, x, z) are coefficients, and lmax is the maximum of the azimuthal quantum number on which the number of terms of the series solution depends.
Figure 2. Flow geometry and coordinate system. The coordinates of a pair of defects are P (0.3H, 0.5H) and Q (0.7H, 0.5H). A defect with positive strength exists at point P and a defect with negative strength exists at point Q. Line segments mean the initial distribution of the director (Equation (18)).
Since the head and tail of a liquid crystal molecule are indistinguishable, we have . From the parity of the spherical harmonic functions, that is,
the expression (?1)l =1 is obtained, which restricts l to even values. We have non-dimensionalized Equation (1) with 1/D, multiplied the resulting equation by the complex conjugate of spherical harmonic functions, , and integrated over the unit sphere to get the ordinary differential equations with respect to Clm (see Appendix 1). The orthonormality of the spherical harmonic functions 
has been used.
To compute the velocity field, we eliminate p from Equation (6) to obtain the vorticity transport equation
where is the vorticity, and N1 (=τxx − τzz) is the first normal stress difference. Stresses N1, τxz, and τzx are explained in the Appendix 2.
We have used the finite-difference scheme to discretize the equations and the implicit Euler method for time integration. Equation (15) is solved using an iterative procedure with the convergence criterion that the average relative error at each node is less than 10−5. We use periodic boundary conditions.
Values of the physical quantities are the fluid density ρ = 103 kg/m3, the absolute temperature T = 320 K, the rotary diffusivity in an isotropic state D = 5.2 × 103 s−1, the number density of molecules c = 2.25 × 1024 m−3, the drag coefficient of rotary molecules ζ = 8.89 × 10−24 kgm2/s, and the side length of the computational domain H = 1 μm. The computational parameters we select are the nematic potential intensity U = 5.0, 5.5, and 6.0, and the long-range order effect (=li/H) = 0.02 − 0.1 with a step of 0.01. The choice of lmax is important, and we have determined it by accounting for both the computational accuracy and load. We report the details in the following chapter.
The mesh size is set to , and the time step Δt* (=ΔtD) is varied depending on ; for example, Δt* = 5 × 10−4 when .
2.3. Initial Values
The initial velocity vector v is 0. The initial values of Clm are obtained as follows: We multiply Equation (12) by , integrate over the unit sphere, and use Equation (14) to get
Assuming that f is in an equilibrium state (no flow) at t* = 0, f is expressed as 
We use the denominator to normalize f in Equation (17). Since f has uniaxial symmetry in the absence of both flow and external field, the order parameter tensor, Equation (4), is rewritten as
where n is a unit vector describing the average local molecular orientation, called the director, and S is the scalar order parameter ranging from 0 in a random orientation state to 1 in a perfect alignment, defined as
The symbol “:” means the double dot product of two tensors. Substituting Equation (3) and Equation (18) into Equation (17), expanding exponential terms into a power series, and expressing the power by the spherical harmonic functions give
where , , 0!! = 1 (n is 0 or even number). In Equation (20), n and S are functions of x and z. We have determined both values (see Appendix 3).
3. Results and Discussion
3.1. Determination of lmax
Before computing the velocity and orientation fields, we must determine the value of lmax. Let us consider a liquid crystal to which simple shear flow is applied. Setting li = 0, Equation (3) reduces to
which is the Maier-Saupe nematic potential  . The coordinate system is shown in Figure 3. The flow is in the x direction, and the velocity gradient is in the z direction. The polar and azimuthal angles of the unit vector u are θ and ϕ, respectively. For simple shear flow, the velocity gradient tensor is expressed as
where is the shear rate. Equation (11), Equation (21), Equation (22) are substituted into Equation (1) to obtain the ODF f. Equation (20) can be used as an initial value of Clm, but it is a function only of time t. At t* = 0 the director n is assumed to orient in the x direction, so that we set θm = π/2 and ϕm = 0, where θm and ϕm are the polar and azimuthal angles of the director, respectively.
The azimuthal angle of the director ϕm is always 0 because the director is in the x-z plane in Figure 3. The polar angle θm is obtained from S as follows:
Figure 4 shows the transient behaviors of θm and S at lmax = 4 − 10 (no convergent solution was obtained at lmax = 2) when the potential intensity U = 5 and the dimensionless shear rate and 40. In Figure 4(a), the director rotates and the order parameter oscillates about the equilibrium value S = 0.615. The behaviors at lmax = 8 are almost the same as those at lmax = 10, and the behaviors at lmax = 6 have smaller periods compared with those at lmax = 8 and 10. At lmax = 4, steady-state instead of periodic behaviors are obtained. In Figure 4(b), the director and the order parameter reach their steady state values in a short time at lmax = 6, 8, and 10, whereas the predictions at lmax = 4 are even qualitatively different from the results at other lmax. Thus, lmax = 4 is not available, lmax = 6 is qualitatively acceptable, and lmax ≥ 8 is quantitatively
Figure 3. Coordinate system.
Figure 4. Transient behavior of preferred angle θm and order parameter S for the potential intensity U = 5. No convergent solution was obtained at lmax = 2. The lines at lmax = 10 overlap with those at lmax = 8. (a) . The director rotates and the order parameter oscillates about the equilibrium value S = 0.615; (b) . Both the director and the order parameter reach their steady state values.
satisfactory. Expansion terms required to approximate the ODF depend on the potential intensity and shear rate. When U and are large, more terms are necessary because the ODF becomes steep and non-uniaxial.
3.2. Orientation Fields
Since the flows induced by annihilation of a pair of defects are not large, and the selected potential intensity U in the present study is not high (5, 5.5, 6), we have used lmax = 6 for the following computations.
Figure 5 shows a sequence of predictions of transmitted light intensity observed under crossed nicols for . The transmitted light intensity I is evaluated using
Here, I0 is the incident light intensity, and θm is the polar angle of the director obtained by Equation (23). A pair of defect cores attracts each other with time, and brushes connecting the defects become short. Finally the defects are annihilated at t* = 3.85. Since the initial orientation profile is symmetric with respect to z* = 0.5, the defect cores move along this line. After finishing the annihilation, the orientation state does not become homogeneous immediately, but slightly gray areas are discernible at t* = 5. A completely homogeneous orientation state is achieved at t* = 6.05.
Figure 6 shows the profiles of S for the same parameters and times as
Figure 5. A time series of the molecular orientation fields at . The transmitted light intensity is numerically predicted. The horizontal axis is x* (=x/H) and the vertical axis is z* (=z/H). The polarizing and analyzing directions are ±45˚ with respect to the x* axis. A field where the director orients to the polarizing or analyzing directions is dark, and a field where it orients parallel to the x* or z* axis is bright, while the other field is gray.
Figure 5. Since the order parameters near the defect cores are low, the profile of S looks like an inverted cone. The two inverted cones gradually approach and coalesce laterally with each other as time advances. The two sharp ends merge into a single end at t* = 3.85, which coincides with the time the defects
Figure 6. Time series of the order parameter fields at .
annihilate in Figure 5. Just after the annihilation, we can see slight inhomogeneity in S.
Figure 7(a) plots the locations of a pair of defect cores before annihilation for several . The defect cores approach at a constant speed, but they slow down just before the annihilation. The defect core with the strength s = +1/2 (at point P in Figure 2) moves slightly faster than that with the strength s = −1/2 (at point Q in Figure 2). Thus, the annihilation takes place at a position slightly over x* = 0.5.
Figure 7(b) shows the relationship between the time required for the
Figure 7. (a) Time series of locations of a pair of defects. Symbols “×” are the positions at which the two defect cores coalesce with each other; (b) Relationship between the time required for the annihilation process and .
annihilation process and the long-range order parameter . A large gives a short annihilation time because the attractive force between defects becomes strong for large . For example, at , and at (refer to Figure 5 and Figure 6); that is, increasing by a factor of 5 decreases by a factor of 50. Therefore, has a large effect on annihilation.
3.3. Velocity Fields
Figure 8 shows the velocity distributions at t* = 1, 3.85, and 5 for as well as the locations of defect cores. An arrow at the bottom of each figure is the reference velocity vector. It is confirmed that flows are induced during annihilation of a pair of defects. Since the molecular orientation state is symmetric with respect to z* = 0.5, the velocity distribution is also symmetric. Flows to the right are observed in the vicinity of z* = 0.5. As a result, a counterclockwise vortex in the upper half region and a clockwise one in the lower half region are generated. Four vortices (two in each half region) are generated at t* = 1 and 3.85. By comparing the velocity vectors near the left and right defect cores at t* = 1, we find that the velocity vector and vortex on the left side are larger, and that the velocity distributions depend on the defect structure. Both right and left vortices rotate in the same direction at t* = 1, so that they coalesce into a larger vortex at t* = 3.85. Since the induced velocity vector is maximum spatially on the line connecting the defect cores (z* = 0.5), it is efficient to use the flow between defects when using the flow induced during their annihilation. We have checked the velocity profiles at times different from those of Figure 8 to verify that the velocity vector has its maximum value when annihilation finishes (t* = 3.85).
We define vmax as the maximum velocity value in the computational region at the time when annihilation finishes. Figure 9(a) plots the relationship between vmax and . When increases, the induced velocity also increases. However, the effect of on the velocity is not large compared with that on the annihilation
Figure 8. Velocity profiles for . Red circles are the locations of defect cores.
Figure 9. Relationship between maximum velocity vmax and and U. (a) Effect of ; (b) Effect of U for
time ; the relationship between vmax and is almost linear. Figure 9(b) plots vmax against the nematic potential intensity U for . The induced velocity increases with increasing U. We explain this result as follows: Since annihilation stems mainly from the spatial gradient of S, an increase in U corresponds to an increase in S, resulting in a steep spatial gradient for S. Thus, a liquid crystal with higher liquid crystallinity (a low-temperature liquid crystal) may generate a faster flow.
In this study we have predicted flows induced by the annihilation of a pair of liquid crystal defects using the Doi theory with the Marrucci-Greco potential and the constitutive equation of Feng et al. The long-range order effect on the time required for the annihilation process of a pair of defects is remarkable; when the long-range order is large, the annihilation time becomes short. We have shown that a flow is induced by the annihilation and that several vortices are generated in the vicinity of the defects. The maximum flow is obtained on the line connecting the two defect cores in space, and at the time, the annihilation is just finished. The maximum value of the induced velocity is on the order of 10 μm/s in our study. The induced velocity becomes large when the long-range order and nematic potential strength are high.
This work was partially supported by the Japan Society for the Promotion of Science KAKENHI (Grant No. 25289035).
The time evolution equation for coefficients Clm is expressed as
The non-zero components of the order parameter tensor S and fourth-order tensor Q are expressed in terms of Clm as
The first normal stress difference N1 and the shear stresses τxz and τzx are expressed in terms of S and Q as
The difference between τxz and τzx is only the underlined term.
With respect to the director n, let us define the angle between n and the x axis as θc, so that the angular momentum equation of the Leslie-Ericksen theory with the one-constant approximation of the elastic constants in the molecular field reduces simply to
in an equilibrium state  . When defects are included, a fundamental solution of Equation (A1) is
where s is the defect strength. Since Equation (A1) is linear, the superposition of solutions is effective. Thus, when a defect with s = +1/2 exists at P (x1, z1) and a defect with s = −1/2 exists at Q (x2, z2) in Figure 2, the director distribution around the defect cores is
Finally, we have modified the values of Equation (A3) so that they fit the periodic boundary condition. We denote the initial distribution of the director by line segments in Figure 2.
For the scalar order parameter, we set S = 0 at defect cores and S = Seq at the other region. Seq depends only on the nematic potential intensity U, and is obtained from Equation (1) without flow terms. For example, Seq = 0.615 at U = 5.