Order Reconstruction in a Nanoconfined Nematic Liquid Crystal between Two Coaxial Cylinders

The dynamics of a disclination loop (s = ±1/2) in nematic liquid crystals constrained between two coaxial cylinders were investigated based on two-dimensional Landau–de Gennes tensorial formalism by using a finite-difference iterative method. The effect of thickness (d = R2 − R1, where R1 and R2 represent the internal and external radii of the cylindrical cavity, respectively) on the director distribution of the defect was simulated using different R1 values. The results show that the order reconstruction occurs at a critical value of dc, which decreases with increasing inner ratio R1. The loop also shrinks, and the defect center deviates from the middle of the system, which is a non-planar structure. The deviation decreases with decreasing d or increasing R1, implying that the system tends to be a planar cell. Two models were then established to analyze the combined effect of non-planar geometry and electric field. The common action of these parameters facilitates order reconstruction, whereas their opposite action complicates the process.


Introduction
The equilibrium configuration of a confined nematic liquid crystal (NLC) is caused by interplay among surface interaction, elastic distortion, and finite-size effect. Confined systems that suffer from continuous symmetry-breaking transitions often display topological defects [1][2][3][4], which are utilized in new-generation LC devices. Therefore, study of defects is an important field in physics. In particular, defects in nanoconfined LC cells have received increased research attention. A weak local perturbation in an LC confined in a nanoscale cell can stimulate an apparent mesoscopic or even macroscopic response [5].
Eigenvalue exchange/order reconstruction occurs under severe confinement in the presence of antagonistic anchorings, and this mechanism relaxes surface-induced frustrations. Order reconstruction structures have been investigated in detail [6][7][8][9][10][11] for various boundary conditions in nematic cells bounded by parallel walls, for which the characteristic linear size of the confining plates is large and virtually infinite compared with the cell thickness. These structures are stable in a sufficiently thin cell [6,7], whose gap is comparable with the order parameter length ξ 0 (characteristic length for order-parameter changes), or when a sufficiently strong electric or magnetic field is applied [8,9]. Studies [12][13][14] revealed that mesoscopic Landau-de Gennes theory can accurately predict the structural and phase behavior of severely confined LCs.
Although LC defects have been intensively studied theoretically and experimentally, LC shells have been rarely investigated. LC shells have received increased research interest because of their potential to generate colloids with a valence, which can be used to build colloidal architectures for photonic applications. Several scholars have numerically simulated the possible defect configurations in nematic and smectic shells [15][16][17].
NLCs in a cylindrical cavity have been studied for approximately 40 years. Cylindrical geometry is the most convenient geometry after planar cell geometry for theoretical and experimental analyses. As such, all experimentally observed and theoretically possible configurations of the director field Ñ n yield a complete or approximate analytic description [18]. The possible equilibrium configurations of an NLC in a cylindrical cavity primarily depend on how Ñ n is anchored to the cylinder surface. In this study, we investigated the structural behavior of a nanoscale NLC system confined in a cylindrical cavity with a topological defect loop by using the finite-difference iterative method. The outline of the paper is as follows. In Section 2, we introduce the phenomenological model employed and describe the geometry of the problem and our parametrization. The results are presented in Section 3, and the conclusions are summarized in Section 4.

Free Energy
Our theoretical argument is based on Landau-de Gennes theory [19], wherein the orientational order of LC is described by a second-rank symmetric and traceless tensor [3]: where Ñ e i is the orthogonal unit vector representing the eigenvector of Q, and λ i is its eigenvalue.
The range of eigenvalues must be λ i Pˆ´1 3 , 2 3˙t o interpret Q as the traceless second-moment tensor of the molecular distribution function. Q vanishes in the isotropic phase, whereas this parameter has two degenerate eigenvalues in the uniaxial ordering and can be represented by where Ñ n is the nematic director pointing along the local uniaxial ordering direction, and S is the uniaxial scalar parameter expressing the magnitude of fluctuations about the nematic director. In Equation (2), S can be either positive or negative. The ensemble of molecules represented by Q tends to align along Ñ n when S is positive and tends to lie in the plane orthogonal to Ñ n when S is negative. Finally, the LC is in a biaxial state when all the eigenvalues of Q are distinct. The degree of biaxiality is expressed by the parameter β 2 defined as [20] This equation is a convenient parameter for illustrating spatial inhomogeneities of Q and ranges in the interval [0,1]. All uniaxial states with two degenerate eigenvalues correspond to β 2 = 0, whereas states with maximal biaxiality correspond to β 2 = 1. As tr(Q 3 ) = 3detQ, states with β 2 = 1 are those with det Q = 0, which implies that at least one eigenvalue of Q vanishes.
The Landau-de Gennes free energy density of LC is given by Hence, Ñ e ϑ is always an eigenvector of Q. This configuration rules out distortions twisted along Ñ e ϑ . Q rz in Equation (9) reveals the departures of the eigenframe´Ñ e 1 , Ñ e z¯, with both frames coinciding when Q rz = 0. We denote the free boundary conditions at the upper and lower lateral walls z =˘d z /2, with initial uniaxial ordering by total rotations of´π/2 (z = d z /2) and +π/2 (z =´d z /2). These boundary conditions are consistent with the generation of a defect loop with topological charge s =´1/2. Figure 1b shows the director profile in a cross-section along an arbitrary radius of the system. The results are also applicable to the defect loop with topological charge s = 1/2, with exchanged conditions of the upper and lower walls.

4
( ) sin ,0,cos n = β β  , where â is the angle with respect to z e  (as shown in Figure 1) and the order tensor Q can be represented as Hence, e ϑ  is always an eigenvector of Q. This configuration rules out distortions twisted along e ϑ  . Qrz in Equation (9)  We denote the free boundary conditions at the upper and lower lateral walls z = ±dz/2, with initial uniaxial ordering by total rotations of −π/2 (z = dz/2) and +π/2 (z = −dz/2). These boundary conditions are consistent with the generation of a defect loop with topological charge s = −1/2. Figure 1b shows the director profile in a cross-section along an arbitrary radius of the system. The results are also applicable to the defect loop with topological charge s = 1/2, with exchanged conditions of the upper and lower walls.

Scaling and Dimensionless Evolution Equations
We introduced the following dimensionless quantities: is the superheating order parameter at the nematic superheating temperature T**, and is the characteristic length for order-parameter changes. Equations (4), (5

Scaling and Dimensionless Evolution Equations
We introduced the following dimensionless quantities: B{4C is the superheating order parameter at the nematic superheating temperature T**, and ξ 0 " is the characteristic length for order-parameter changes. Equations (4), (5), and (8) can be reduced to where the reduced parameter r A " 24AC B 2 defines the temperature scale. Isotropic-nematic transition occurs at r A " 8{9. In polar cylindrical coordinates, Equation (11) can be expressed precisely by Thus, the reduced uniaxial ordering has the form where r s " S eq q 0 " 1`a1´r A is the reduced uniaxial scalar parameter at equilibrium.
We compute the evolution of LC with dynamic theory for tensor order-parameter field Q(r,z,t). The local values of the scalar-order parameter S and the director Ñ n can be calculated from Q by using the highest eigenvalue and the associated eigenvector, respectively. According to [23], the evolution equation describing the dynamics of Q can be written as where Γ = 6D*/[1´3tr(Q 2 )] 2 , D* is the rotational diffusion for the nematic, and δ f {δQ is assumed to be symmetrical. Numerical calculations are performed using the reduced variables. When the functional derivatives in Equation (15) are evaluated, the following evolution equations for r Q can be obtained.
r Q rr`p r Q 2 q rr´´t r r Q 2¯˜r Q rr 4`1 3¸´2 r r 2´r Q rr´r Q ϑϑr Q rr,r rr r`1 r r r Q rr,r r`r Q rr,r zr z  , with r Γ " ΓˆpBq 0 q. We adopt the two-dimensional finite-difference method developed in our previous studies [24,25] to obtain the numerical simulation results. Here, we let the system relax from the initial boundary conditions given in Section 2.2, and the initial conditions of the upper and lower half in the bulk are specified by total rotations of´π/2 and +π/2, respectively, consistent with the boundary conditions at the upper and lower walls given above. We consider only the equilibrium configurations that correspond to the global minimum F. In our numerical calculations, a discretization with a time step given by 10´1 0 is sufficient to guarantee the stability of the numerical procedure. In addition, our equilibration runs are 2ˆ10 6 , which is adequate for the system to reach the equilibrium.

Results
According to parameters given in [26], a = 0.195ˆ10 6 J/m 3 , B = 7.155ˆ10 6 J/m 3 , C = 8.82ˆ10 6 J/m 3 , L = 10.125ˆ10´1 2 J/m can be easily obtained. In our simulations, the scaled temperature is set to r A " 2{3, which corresponds to r s " 1`a1{3. The rotational diffusion D* is set to 0.35, which is the value used in [26], and the dielectric anisotropy is set to ∆ε = 13.2 (ε K = 6.5). We then obtain the exact values ξ 0 = 2.64 nm and E 0 = 43.1 V/um. We focus on the dependence of structure transition on thickness (d) and then on the structure transition induced by electric field ( Ñ E ).

Structure Transition with Different Values of Thickness d and Internal Radius R 1
We first analyze the equilibrium textures in the case of R 1 = 30 ξ 0 with different values of d. Figures 2 and 3 show the director field profile and the calculated biaxiality β 2 in a cross-section along an arbitrary azimuth for different values of d. At relatively large thicknesses, the nematic system shows the defect core structure (Figure 2a,b) and is uniaxial everywhere (β 2 = 0), except for a small region around the defect core (Figure 3a,b). When d decreases, the defect core clearly explodes along z (Figure 2c), and a large biaxiality propagates along z inside the system (Figure 3c). Further decrease in d results in the creation of a biaxial wall, which connects the two orthogonal uniaxial directions imposed by two cylindrical surfaces (Figures 2d and 3d) at a critical value of d c = 11.2 ξ 0 . In summary, the system has transitioned from the eigenvector rotation configuration with a defect loop into the eigenvalue exchange/order reconstruction configuration by developing a thin biaxial nematic layer, which forms a cylindrical wall in the nematic system. The transition process is similar to that conducted in a planar cell in our previous studies [24,25]. Figures 2 and 3 show that the defect center is not located in the middle of the system but shifts to the internal surface; this phenomenon differs from the case of a planar cell [24,25]. Non-planar cylindrical geometry causes disclination loops to shrink, and the loop will not disappear under internal surface confinement [27]. Figure 3 shows that the deviations ∆(d) for different thicknesses (d) are 1.6 ξ 0 (d = 15 ξ 0 ), 1.3 ξ 0 (d = 13 ξ 0 ), 0.95 ξ 0 (d = 11.8 ξ 0 ), 0.85 ξ 0 (d = 11.2 ξ 0 ), 0.45 ξ 0 (d = 9 ξ 0 ), and 0.18 ξ 0 (d = 6 ξ 0 ). These finding indicate that ∆(d) decreases with decreasing d, even after eigenvalue exchange (Figure 3e,f). This behavior is due to the gradual transformation of the system to a planar structure with decreasing d.
The relative deviations δ(d) = ∆(d)/d are calculated and normalized by δ(d)/δ(15 ξ 0 ) to analyze the influence of non-planar geometry. Figure 4 shows the curves of δ(d)/δ(15 ξ 0 ) as a function of d/ ξ 0 . In the figure, the relative deviations decrease with decreasing thickness d.
A clearer explanation of this phenomenon is as follows: For a planar hybrid-aligned cell (R 1 Ñ8), the defect center is located in the middle of the system, with the energy on both sides of the defect equal. While for the non-planar cylindrical geometry, if the director configuration keeps the configuration of a planar cell, the energy on the inner side of the defect will be lower than that on the outer side, because that the volume on the inner side is smaller, relative to the outer side. Thus the defect center shifts to the inner side, and the director configuration of the nematic meanwhile changes, until the system reaches equilibrium. For fixed R 1 , the volume differences on the two sides of the middle of the system decrease with decreasing d, thus deviations decrease with decreasing d.      The energy calculation shows that the energy in the "inner" region of the cylindrical geometry is higher than that in the "outer" region, indicating that more energy can be saved by shrinking the disclination to a smaller radius to reduce its total length. In addition, the shrinking of the disclination to a smaller radius also helps to decrease the total energy of the system.  The energy calculation shows that the energy in the "inner" region of the cylindrical geometry is higher than that in the "outer" region, indicating that more energy can be saved by shrinking the disclination to a smaller radius to reduce its total length. In addition, the shrinking of the disclination to a smaller radius also helps to decrease the total energy of the system. The energy calculation shows that the energy in the "inner" region of the cylindrical geometry is higher than that in the "outer" region, indicating that more energy can be saved by shrinking the disclination to a smaller radius to reduce its total length. In addition, the shrinking of the disclination to a smaller radius also helps to decrease the total energy of the system.  The conditions of R 1 = 20 ξ 0 and R 1 = 40 ξ 0 are also investigated to explore the influence of internal radius R 1 on equilibrium texture and transition. The critical d c values for different R 1 values are shown in Figure 5. The value of d c increases with decreasing R 1 . This behavior indicates that the smaller the internal radius R 1 , the earlier the biaxial transition occurs. Therefore, a non-planar cylindrical system induces easy biaxial transition. The cylindrical system gradually flattens with increasing R 1 . Thus, we predict that d c gradually becomes equal to that of a planar cell with increasing R 1 . Figure 6 shows the biaxiality β 2 in a cross-section for different values of R 1 with fixed thickness d = 11.8 ξ 0 . Deviations ∆(d) are 1.25 ξ 0 (R 1 = 20 ξ 0 ), 0.95 ξ 0 (R 1 = 30 ξ 0 ), and 0.8 ξ 0 (R 1 = 40 ξ 0 ); this finding indicates that ∆(d) decreases with increasing R 1 or with gradual flattening of the system. That is because for fixed d, the volume ratio on the two sides of the middle of the system decreases with increasing R 1 . We predict that ∆(d) gradually approaches zero when R 1 approaches infinity. The conditions of R1 = 20 ξ0 and R1 = 40 ξ0 are also investigated to explore the influence of internal radius R1 on equilibrium texture and transition. The critical dc values for different R1 values are shown in Figure 5. The value of dc increases with decreasing R1. This behavior indicates that the smaller the internal radius R1, the earlier the biaxial transition occurs. Therefore, a non-planar cylindrical system induces easy biaxial transition. The cylindrical system gradually flattens with increasing R1. Thus, we predict that dc gradually becomes equal to that of a planar cell with increasing R1. Figure 6 shows the biaxiality β 2 in a cross-section for different values of R1 with fixed thickness d = 11.8 ξ0. Deviations Δ(d) are 1.25 ξ0 (R1 = 20 ξ0), 0.95 ξ0 (R1 = 30 ξ0), and 0.8 ξ0 (R1 = 40 ξ0); this finding indicates that Δ(d) decreases with increasing R1 or with gradual flattening of the system. That is because for fixed d, the volume ratio on the two sides of the middle of the system decreases with increasing R1.We predict that Δ(d) gradually approaches zero when R1 approaches infinity.   The conditions of R1 = 20 ξ0 and R1 = 40 ξ0 are also investigated to explore the influence of internal radius R1 on equilibrium texture and transition. The critical dc values for different R1 values are shown in Figure 5. The value of dc increases with decreasing R1. This behavior indicates that the smaller the internal radius R1, the earlier the biaxial transition occurs. Therefore, a non-planar cylindrical system induces easy biaxial transition. The cylindrical system gradually flattens with increasing R1. Thus, we predict that dc gradually becomes equal to that of a planar cell with increasing R1. Figure 6 shows the biaxiality β 2 in a cross-section for different values of R1 with fixed thickness d = 11.8 ξ0. Deviations Δ(d) are 1.25 ξ0 (R1 = 20 ξ0), 0.95 ξ0 (R1 = 30 ξ0), and 0.8 ξ0 (R1 = 40 ξ0); this finding indicates that Δ(d) decreases with increasing R1 or with gradual flattening of the system. That is because for fixed d, the volume ratio on the two sides of the middle of the system decreases with increasing R1.We predict that Δ(d) gradually approaches zero when R1 approaches infinity.

Structure Transition Induced by Electric Field E 
In this section, two models are established to analyze the combined effect of geometry and electric field E  . Boundaries on both coaxial cylindrical surfaces are exchanged for the two models, as shown in Figure 7. Simulations are conducted at internal radius R1 = 30 ξ0 and thickness d = 15 ξ0.  Figures 8 and 9 show the director field profile and the calculated biaxiality β 2 in a cross-section along an arbitrary azimuth as a function of E  for model I, and Figures 10 and 11 show the results for model II, respectively. The disclination loops shift to the internal surface at 0 E =  , which is caused by the non-planar cylindrical geometry. The defect is pushed further toward the internal surface for model I with increasing E  , whereas the defect initially shifts to the middle of the system and then to the external surface for model II. Because of the strong anchoring boundaries, the defects exhibit dramatic changes in shape as the distance between defect center and the surface boundary decreasing, especially when the defect center lies very close to the surface, a biaxial layer is established, that is order reconstruction occurs [9,24,25] Non-planar cylindrical geometry induces disclination loops to shrink for both models, whereas electric field E  , which tends to enforce the nematic director along the r e  axis, expulses the defect to the internal and external surfaces for models I and II, respectively. It means that the effects of the cylindrical geometry and electric field E  are common for model I, but opposite for model II. Hence, the transition process differs between the two models. The common action of the cylindrical geometry and electric field E  Note that the radial direction uniform electric field has been assumed, which is difficult to implement or almost impossible in practice. Our results only give the qualitative analysis about the effect of the electric field as well as the combined effect of geometry and the electric field. We can  Figures 8 and 9 show the director field profile and the calculated biaxiality β 2 in a cross-section along an arbitrary azimuth as a function of r E for model I, and Figures 10 and 11 show the results for model II, respectively. The disclination loops shift to the internal surface at r E " 0, which is caused by the non-planar cylindrical geometry. The defect is pushed further toward the internal surface for model I with increasing r E, whereas the defect initially shifts to the middle of the system and then to the external surface for model II. Because of the strong anchoring boundaries, the defects exhibit dramatic changes in shape as the distance between defect center and the surface boundary decreasing, especially when the defect center lies very close to the surface, a biaxial layer is established, that is order reconstruction occurs [9,24,25]. Our results show that order reconstruction occurs at critical values of r E c pIq " 0.21 and r E c pIIq " 0.25 for models I and II, respectively. For comparison, a planar cell is also simulated with the same parameters, and the corresponding reduced critical value of electric field is r E c0 " 0.24. The values are clearly ranked as follows: r E c pIq ă r E c0 ă r E c pIIq. Non-planar cylindrical geometry induces disclination loops to shrink for both models, whereas electric field Ñ E , which tends to enforce the nematic director along the Ñ e r axis, expulses the defect to the internal and external surfaces for models I and II, respectively. It means that the effects of the cylindrical geometry and electric field Ñ E are common for model I, but opposite for model II. Hence, the transition process differs between the two models. The common action of the cylindrical geometry and electric field Ñ E for model I makes the defect close to the surface boundary more easily than model II, since the opposite action of the cylindrical geometry and electric field Ñ E makes it difficult for the defect to be close to the surface boundary. From the above, the common action of the cylindrical geometry and electric field Ñ E for model I facilitates order reconstruction, while the opposite action of the cylindrical geometry and electric field Ñ E for model II complicates order reconstruction. For a planar cell, only the action of electric field Ñ E exists. Considering these factors, our result r E c pIq ă r E c0 ă r E c pIIq is reasonable. Note that the radial direction uniform electric field has been assumed, which is difficult to implement or almost impossible in practice. Our results only give the qualitative analysis about the effect of the electric field as well as the combined effect of geometry and the electric field. We can conclude that for the actual cylindrical symmetric radial direction electric field E prq " E 0 r , the same conclusion is reached in nature. The detailed and definite results for the actual cylindrical symmetric radial direction electric field are a task for the future.
, the same conclusion is reached in nature. The detailed and definite results for the actual cylindrical symmetric radial direction electric field are a task for the future.   conclude that for the actual cylindrical symmetric radial direction electric field ( ) 0 E E r r = , the same conclusion is reached in nature. The detailed and definite results for the actual cylindrical symmetric radial direction electric field are a task for the future.

Conclusions
Structural transitions of an NLC confined between two coaxial cylinders were investigated through Landau-de Gennes theory by using a two-dimensional finite-difference iterative method. We considered the cylindrical symmetry configuration about the cylindrical axis. The effects of the shell thickness d, as well as the strength of the electric field E  , were also analyzed. The results show that order reconstruction occurs in a thin enough shell or under a strong enough electric field. The non-planar system causes the defect center to deviate from the middle of the system to the because of the generated defects. The evident effect of Ñ E on the dynamic behavior of LC around the defects and transitions between the topologically different states can contribute to the design of advanced LC modes.
Our simulation describes an isolated and stabilized defect reasonably well. A natural question is whether our simulation describes a realistic scenario within such a coaxial cylindrical nematic system. In fact, defect rings with alternating winding signs exist when the defect spacing in the z direction is large enough. With defect spacing decreasing, the defect rings gradually annihilate at a characteristic length. Identifying the detailed effects of the defect spacing in the z direction as well as the electric field on the defect structure are tasks for the future.
Comparison of experimental works with numerical results in nematic shells is rare because shells are difficult to produce in a controlled manner. Thus far, microfluidics has provided a natural method to overcome this limitation. Nevertheless, future theoretical and experimental investigations must be performed to determine the effects of electromagnetic field, temperature, or boundary anchoring on nematic shells and even smectic or cholesteric shells, especially when degenerate planar anchoring condition is prescribed on the cylinder surface; this would be a good and interesting topic. In addition, the geometry of a bookshelf on curved surfaces is also an interesting task for the future.