Dynamics of Order Reconstruction in a Nanoconfined Nematic Liquid Crystal with a Topological Defect

At the wall in a hybrid nematic cell with strong anchoring, the nematic director is parallel to one wall and perpendicular to the other. Within the Landau-de Gennes theory, we have investigated the dynamics of s = ±1/2 wedge disclinations in such a cell, using the two-dimensional finite-difference iterative method. Our results show that with the cell gap decreasing, the core of the defect explodes, and the biaxiality propagates inside the cell. At a critical value of dc* ≈ 9ξ (where ξ is the characteristic length for order-parameter changes), the exchange solution is stable, while the defect core solution becomes metastable. Comparing to the case with no initial disclination, the value at which the exchange solution becomes stable increases relatively. At a critical separation of dc ≈ 6ξ, the system undergoes a structural transition, and the defect core merges into a biaxial layer with large biaxiality. For weak anchoring boundary conditions, a similar structural transition takes place at a relative lower critical value. Because of the weakened frustration, the asymmetric boundary conditions repel the defect to the weak anchoring boundary and have a relatively lower critical value of da, where the shape of the defect deforms. Further, the response time between two very close cell gaps is about tens of microseconds, and the response becomes slower as the defect explodes.


Introduction
Topological defects arise as a result of broken continuous symmetry and are ubiquitous in nature, from microscopic condensed matter systems governed by quantum mechanics to a universe in which gravity plays a decisive role [1][2][3]. Defects in liquid crystals (LCs) have been the subject of much interest, still offering unsolved problems. Commonly observed defects in the uniaxial nematic phase are typically point defects with topological charge s = 1 and line defects with topological charge s = ±1/2 [4]. There are two types of s = ±1/2 disclination lines [5]. The first type is wedge disclination, with the rotation vector parallel to the disclination line, while the second is twist disclination, with the rotation vector perpendicular to the disclination line. The region where the presence of a defect causes apparent deviations from bulk ordering is referred to as the defect core [6].
Eigenvalue exchange (also called order reconstruction) was first shown by Schopohl and Sluckin within the core of s = ±1/2 wedge disclinations [7] and has aroused numerous subsequent studies on the detailed core structure in nematic LCs for various boundary conditions theoretically [8][9][10][11][12][13] and experimentally [14][15][16][17]. In particular, Palffy-Muhoray [8] investigated the eigenvalue exchange problem in a hybrid alignment nematic (HAN) cell with perpendicular alignment of NLC molecules on one surface and parallel alignment on the other. A structure transition from the eigenvector rotation configuration to the eigenvalue exchange is expected when the cell gap is decreased. Similar results are achieved in a twisted nematic (TN) cell with π/2 boundary twist configuration [9].
Recently, the atomic force microscope (AFM) has been used as a unique tool to measure the separation dependence of surface forces with unprecedented accuracy and flexibility [18][19][20][21][22][23][24][25][26]. This is not only important for studying the molecular and surface forces due to their potential application in self-assembly at the nanoscale, but is also important for the structural analysis of the interfaces of soft matter and complex fluids in contact with a solid wall. As an example, the structure of LC layers have been successfully analyzed using the force spectroscopy mode of an AFM [18][19][20]. Particularly, Carbone et al. [18] give the experimental evidence of the eigenvalue exchange in a thin HAN layer in the presence of a topological defect. We know that the speed of the AFM probe is about 1-5 µm/s, which is much slower. corresponding to the dynamic response of the structural transition (~10-100 µs) [14][15][16], in other words, the speed of the AFM probe is slow enough for the system to undergo a structural transition; so, the AFM can be used to investigate the dynamic process of the structural transition as the call gap varies.
Recently, many theoretical investigations have focused on the dynamical behavior of structure transition, using the Q tensor numerical model [27][28][29][30][31], which also describe the effect of order reconstruction-eigenvalue exchange. In this study, we will investigate the dynamics of biaxial transition in a nanoconfined HAN layer with a topological defect numerically, using the two-dimensional finite-difference iterative method. Our study is based on the Landau-de Gennes theory [32] describing the orientational order of an LC in terms of a second-rank tensor Q, which encompasses both uniaxial and biaxial states.

Free Energy
Our theoretical argument is based on the Landau-de Gennes theory [32], in which the orientational order of an LC is described by a second-rank symmetric and traceless tensor [33]: (1) Here, λ i and e i are the i th eigenvalue and the i th eigenvector of Q, respectively. In the isotropic phase, Q vanishes. In the uniaxial ordering, Q has two degenerate eigenvalues and can be represented by (2) where S is the uniaxial scalar parameter, and the unit vector, n  , is the nematic director pointing along the local uniaxial ordering direction. In Equation (2), S can have either sign: when it is positive, the ensemble of molecules represented by Q tends to be aligned along n  , whereas when S is negative, it tends to lie in the plane orthogonal to n  .
Finally, when all eigenvalues of Q are distinct, the LC is in a biaxial state. The degree of biaxiality is measured by the biaxiality parameter, β 2 , defined as [34] (3) which is a convenient parameter for illustrating spatial inhomogeneities of Q and ranges in the interval [0, 1]. In all uniaxial states with two degenerate eigenvalues, β 2 = 0, while states with maximal biaxiality correspond to β 2 = 1. Since tr(Q 3 ) = 3detQ, the states with β 2 = 1 are precisely those where detQ = 0, which further implies that at least one eigenvalue of Q vanishes.
Following the notation in [7], the Landau-de Gennes free energy density of LC is given by is the bulk energy that describes a homogeneous phase. In f bulk , B and C are positive constants, and A is assumed to vary with temperature T in the form A = a (T − T*), where a is a positive constants and T* is the nematic supercooling temperature. Equation (4) gives the bulk equilibrium value of the uniaxial scalar order parameter in Equation (2), = − (1 + 1 − ) , which also depends on the temperature. The free-energy, f e , which penalizes gradients in the tensor order parameter field, is given in the form The coefficients, L i , are related to the splay, twist and bend elastic constants. In this work, we use the one-elastic-constant approximation just for simplicity, where the splay, twist and bend elastic constants have a common value that depends quadratically on the scalar order parameter K = 2S 2 L 1 , and Equation (5) reduces to (6) The boundary conditions at the LC interfaces are taken into account by the surface term, f s , describing the interaction between the nematic molecules close to the substrate and the substrate itself, which is given by (7) where W s = W/S is the anchoring strength, with W being the anchoring strength in the Frank elastic theory, and Q s is the value of the tensor order parameter preferred by the surface [27].
The stable solution of Q on the boundary condition must satisfy the Neumann condition (8) where ν k is the k th component of the outward vector normal to the substrate. This expression allows us to impose on the surface substrate uniaxial, as well as biaxial conditions. For rigid anchoring, this is equivalent to imposing the Dirichlet condition on the surface, i.e., After an appropriate rescaling of the variables, we can simplify the calculating process. Here, we follow the rescaling of Schopohl and Sluckin [7] by defining the following dimensionless quantities: ≡ /[ (9 ) ] , ≡ , ≡ , where = − is the order parameter at the isotropic-nematic phase transition point and = = is the characteristic length for order-parameter changes; then, the rescaled form of Equations (4)- (7) can be written as where the rescaled parameter, = , defines the temperature scale and = , with = − .
We notice that the isotropic-nematic transition takes place at = .
In the scaling employed above, the rescaled uniaxial ordering has the form is the rescaled uniaxial scalar parameter at equilibrium.

Geometry of the Problem
We choose a HAN cell with strong anchoring boundary conditions; the nematic director is parallel to one wall and perpendicular to the other (see Figure 1). The model has two stable degenerate configurations with a rotation difference of π. The plates are placed at z = ±d/2 of a Cartesian coordinate system. The lengths, d x and d y , of the cell along the x-and the y-axes are much We study the structure of s = ±1/2 wedge disclinations. The disclination line is parallel to the y-axis of the simulation cells, and we seek solutions independent of y.
At the two plates z = ±d/2, we enforce the uniaxial anchoring represented by φ −d/2 = 0 and φ +d/2 = π/2, that is (14) (15) On the lateral walls at x = ±d x /2, we also prescribe fixed boundary conditions with uniaxial ordering. The total rotation at x = d x /2 is −π/2, while at x = −d x /2 is π/2. These boundary conditions are compatible with the generation of line defects with topological charge s = −1/2, with which we shall also apply to line defects with topological charge s = 1/2, with the conditions at the lateral walls exchanged.  ) ,

Numerical Methods
We compute the evolution of the LC with a dynamic theory for its tensor order-parameter field, Q(r,t). The local values of the scalar order parameter, S, and the director, n  , can be obtained from Q through its largest eigenvalue and its associated eigenvector, respectively. According to [17], the evolution equation describing the dynamics of Q can be written as (16) The coefficient, Γ, is given by where D* is the rotational diffusion for the nematic. In these equations, it is assumed that δf/δQ has been symmetrized. The stable solution of Q must satisfy the system Equations (16) and, on the boundary, the Neumann condition (Equation (8)) or Dirichlet condition (Equation (9)).
The numerical calculations are to be intended with respect to the scaled variables. When the functional derivatives in Equation (16) are evaluated, one obtains a partial differential equation for in the bulk, with the boundary conditions, where Γ = Γ × (− ) and ∇ = . In general, Equation (17), subject to boundaries Equations (14) and (15), has more than one solution. One is that the structure satisfies the initial conditions given above, i.e., the configuration with a line defect; the second one is a pure HAN configuration consisting of a rotation of the eigenvector of Q (the director) with the z-axis only, i.e., a HAN configuration with no defect; and the third is the configuration with a biaxial layer, where there is little or no net rotation of the eigenvector, i.e., eigenvalue exchange configuration. In this work, we focus on the equilibrium configuration and dynamics with the initial conditions given in Section 2.2.
We use a two-dimensional finite-difference iterative method employed in our previous studies in [35]. To discretize Equation (17), we replace the derivatives with finite difference as follows: (20) In our numerical calculations, we have found that a discretization with a time step given by 4 × 10 −10 is sufficient to guarantee the stability of the numerical procedure. In addition, our equilibration runs take 10 6 , which has been confirmed as sufficient for the system to reach equilibrium state.
We let the system relax from an initial condition under the boundary conditions given in Section 2.2. We have calculated the tensor Q after the system to reach equilibrium state, and then, the current tensor is diagonalized. Accordingly, the director is identified by the eigenvector possessing the largest eigenvalue. Then, the degree of biaxiality is measured by the biaxiality parameter, β 2 , given by Equation (3), and the corresponding free energy is calculated. When the equilibrium state with a different cell gap is reached, we compute the dynamic evolution process between different cell gaps with the evolution equation given by Equation (16).

Results
In this section, we present our numerical results. According to the parameters given in [36] (3)), inside the cell for different values of cell gap d. We find that at relatively large separations, the core structure is weakly influenced by d (see Figure 2a,b), and the nematic is uniaxial everywhere (β 2 = 0), except for a small region around the defect core (Figure 3a,b). Its size is a few characteristic lengths, ξ. When the gap becomes comparable to ξ, the core of the defect explodes obviously along the x direction ( Figure 2c) and a large biaxiality propagates inside of the cell (Figure 3c). The further decrease of the cell gap results in the creation of a biaxial wall, connecting the two orthogonal uniaxial directions imposed by the boundary conditions (Figures 2d and 3d) at a critical value of d c ≈ 6ξ. In a word, the system has transited from the eigenvector rotation configuration with a topological defect into the eigenvalue exchange configuration by developing a thin transitional biaxial nematic layer.
A more detailed analysis of the defect core, directly showing the process of biaxial transition, is given in Figure 4. There, we plot β 2 across the defect center along z = 0 for different values of cell gap d. As the confining gap decreases, the distance at which β 2 reaches its maximum value β 2 = 1 (i.e., the biaxial region) increased, indicating that the defect core becomes increasingly prolate in the x direction. By comparing the two figures in Figure 4, we see that at the cell gap where the defect core explodes, the changes of the biaxial region is faster.    To quantify the influence of the cell gap on the defect structure, we plot in Figure 5 the change in both the width, r x , and height, r z , of the defect core as the cell gap, d, changes. We formally define r x and r z as the distance along the x-and z-axes, respectively, between the defect center and the point where β 2 = 1. It is shown that r z is only slightly affected, while r x increases rapidly as the defect explodes. At a critical value of about d ≈ 6ξ, the value of r x diverges approximately, signaling a nematic order reconstruction. Therefore, the contours representing the biaxiality β 2 become increasingly prolate in the x direction and, finally, into a plane (see Figure 3).
It seems that the transition to the biaxial wall originates from the biaxiality seed concentrated at the −1/2 defect. The biaxial wall would also be created in the absence of a defect [8,9]. According to [8,9], the two orthogonal uniaxial directors prescribed on the plate (HAN cell or TN cell with a π/2 twist) can be connected either through a director rotation or through an order reconstruction. When the cell gap is greater than a critical value d > d c0 , the rotation solution is stable, while the exchange solution is unstable. When d < d c0 , only the exchange solution exists as a stable configuration; no rotation solution exist. In other words, the rotation solutions merge continuously into the exchange solution at a critical value, d c0 . Figure 5. Variation of the defect core size with different thicknesses. The lengths of r x and r z estimate the core size along the x-and z-axes, respectively.
In order to make the mechanism of the structure transition clear, we make a comparison between the two cases (with and without initial disclination). Figure 6 gives the free energy variation inside the cell for three structures (the structure with a defect core, eigenvalue exchange structure and HAN structure with no defect), where the free energy is calculated from numerical simulation results. Here, the HAN structure plays the role of a case of comparison.  Figure 6 shows that when the gap decreases to a critical value of d c * ≈ 9ξ, the exchange solution is stable, while the defect core solution is metastable. Comparing to the case with no initial disclination, we find that for the configuration with a topological defect, the value at which the exchange solution becomes stable increases relatively. That is because the configuration with a topological defect has larger energy, relative to the pure HAN configuration with no defect. With further decreasing of the cell gap until a critical value, d c , a transition of the defect core solution into the exchange solution happens. Below d c , only the exchange solution exists, but no defect core solution. According to the structure given by Figures 2 and 3, we conclude that for the case with initial disclinations, the value at which only the exchange solutions exist is about d c ≈ 6ξ. Figure 6b presents that at a value of d c ~ 7ξ, the two cases (with and without initial disclination) have nearly the same energy, while the simulated transition value with initial disclination is at d c ~ 6ξ. As mentioned in Section 2.2, we have prescribed fixed boundary conditions on the lateral walls, which will have an impact on the structure. In addition, there may be some calculation error. We speculate that for the two cases (with and without initial disclination), the values at which the structural transition takes place (only the exchange solutions exist) are almost the same, i.e., d c = d c0 .

Dynamical Evolution
In order to make it clear how the configuration changes between two states, we make dynamical calculations of the defect core structure between different cell gaps (10ξ → 9ξ, 9ξ → 8ξ, 8ξ → 7ξ, 7ξ → 6ξ). Figure 7 shows the evolution process of the calculated biaxiality β 2 across the defect center along z = 0 as the cell gap changes. It is shown that if the call gap, d, is greater than the critical value of d c , an initial structure with a defect core cannot change to the eigenvalue exchange structure, even at the value where the exchange structure is stable. The response time is about tens of microseconds, and the response becomes slower as the defect explodes (see Figure 7c).
To analyze the influence of the defect core on the dynamic process, we make dynamical calculations of the exchange configurations between different cell gaps. To aid the comparisons, we give the dynamics as the cell gap decreases. Figure 8 shows the evolution process of the calculated biaxiality, β 2 , along x = 0 as the cell gap changes. It is shown that the response time is about several microseconds, and the response becomes faster as the cell gap decreases. Comparing Figure 7 with Figure 8, we find that the response between two exchange structures is much faster than that between two defect core structures.   Figures 7 and 8 show that above the critical value, d c , there are only defect-defect and exchangeexchange transitions, but no defect-exchange or exchange-defect transitions existing. However, the response of an initial defect configuration with d > d c to d ≤ d c (e.g., 7ξ → 6ξ) will undergo a defect-exchange transition (See Figure 7d). We expect that there are cell gap induced hysteresis effects, which are similar to hysteresis effects induced by the voltage in blue-phase liquid crystal [37]. As the cell gap, d, decreases, the defect core becomes increasingly prolate and undergoes a structural transition to the exchange configuration at a critical value of d c ≈ 6ξ. While the cell gap, d, is increasing from d < 6ξ, the system is in the exchange configuration until d ≈ 9ξ, above which the system will change to the defect core configuration under the effect of thermal disturbance.

Weak Anchoring Boundary Conditions
In order to get the influence of the surface anchoring at the boundary plates on the behavior of topological defects, we focus on the Neumann conditions in this section. First, we study the condition of asymmetric boundaries. We impose strong homeotropic anchoring on the top plate and weak planar anchoring on the bottom plate. Figure 9 shows the director field profile inside a cell of thickness d = 10ξ for different values of at the bottom plate. It is shown that as the anchoring strength, , decreases, the defect glides towards the weak anchoring bottom plate and finally coalesces with the surface layer. This is consistent with the results in [38,39]. Next, we fix the anchoring strength at the bottom plate to = 1, to investigate the influence of the cell gap on the topological defects. Figures 10 and 11 show the director field profile, as well as the change in both the width r x and height r z of the defect core inside the cell for different values of cell gap d. Figure 10 shows that the defect core is not located in the middle of the cell, which is due to the asymmetric boundary conditions. As the cell gap, d, decreases, the line defect moves towards the bottom plane. At a critical value d c (w as ) ≈ 5ξ, the defect is transformed into a surface layer. The surface plane defect is then expelled from the cell when the cell gap reaches d c **(w as ) ≈ 2.5ξ. For d c **(w as ) < d < d c (w as ), the surface biaxial layer bridges two uniaxial stares, one homeotropic in the bulk and the other planer on the boundary. The results are consistent with the results in [12,27].   Figure 11. Variation of the defect core size with different thicknesses for the asymmetric boundaries. The lengths of r x and r z estimate the core size along the x-and z-axes, respectively. In the inset, the small region is enlarged. Here, the strong anchoring boundary plays the role of a case of comparison. Figure 11 shows that the change in r z is slight, while r x increases when the cell gap reaches a critical value, d a , where the defect core begins to lose its circular cross-section and stretches along the x direction. By comparing with the strong anchoring conditions (see the inset in Figure 11), we can conclude that the asymmetric boundary condition has a relatively smaller value of d a . That is because the moving of the defect as the cell gap decreases weakens the frustration on the confining surfaces. Figure 11 also gives that the value, d c (w as ), where the defect is transformed into a surface layer, is smaller than that of strong anchoring conditions. Further, we study the condition of symmetric boundaries. We impose the same value of weak anchoring on both the two boundary plates and studied the effect of anchoring strength on the critical value, d c . In our simulation, the anchoring strength is set equal to = 1 at both plates. Figure 12 shows the director field profile inside the cell for different values of cell gap d, which shows that the structural transition takes place at a critical value of d c (w s ) = 3ξ. Compared with the strong anchoring conditions (i.e., Figure 2), we can conclude that weak anchoring boundary conditions decrease the critical value, d c . That is because the weak anchoring boundary weakens the frustration on the confining surfaces.

Conclusions
Within the Landau-de Gennes theory, we carried out a numerical study on the structure of s = −1/2 wedge disclinations and investigated the dynamics of biaxial transition in a nanoconfined HAN layer, using the two-dimensional finite-difference iterative method. Our results show that with the cell gap decreasing, the core of defect explodes and the biaxiality propagates inside the cell. Below a critical thickness, d c , of a HAN layer with a topological defect, an eigenvalue exchange transition occurs, and the defect is squashed flat into a biaxial layer. Because of the presence of the defect core, the value at which the exchange solution becomes stable increases relative to the pure HAN. The dynamic simulations show that before the critical value of d c , an initial structure with a defect core cannot change to the eigenvalue exchange structure, even at the value where the exchange structure is stable. The response time between two very close cell gaps is about tens of microseconds, and the response becomes slower as the defect explodes. For weak anchoring conditions, plates with strong anchoring repel the defects, while plates with weak anchoring allow the defect to escape through the boundary. Because of the weakened frustration, the weak anchoring boundary condition has relatively lower critical values of d a and d c , corresponding to the value where the defect core begins to lose its circular cross-section and where the defect is transformed into a surface layer, respectively.
According to the parameters given, we can easily get the characteristic length ξ ~ 3.96 nm, and the critical separation is ~20 nm, which is consistent with the experimental result in [18]. Furthermore, the response time also agree with the experimental results (~10-100 µs) in [14]. The consistency of numerical calculation and experiment confirms the correctness and rationality of our theory, which further illustrates that the AFM can be used to carry on dynamics research.
Further, the parameters also give that the unit of anchoring strength is W 0 ≈ 5.7 × 10 −4 J/m 2 , which is reasonable. According to b = k/W, corresponds to an extrapolation length of b ≈ 10 −7 m, which is larger than the cell gap of d = 10ξ ~ 3.96 × 10 −8 m, i.e., b/d > 1. Figure 9 shows that at 1 . 0 = s W , a sample of d = 10ξ presents a homeotropic orientation, which is consistent with the results of [39]. It can be predicted that the critical separation and the response time should be affected by the temperature; the detailed results, as well as the force curve with the Q tensor numerical model are tasks for the future.