Structures and phase transition of liquid crystals in a dynamic slit confinement

Liquid crystals (LCs) have received considerable attention due to their wide range of application in displays, biosensors, optical switches, and many other electro-optical devices. These applications of LCs are based on controlling the orientational and translational order in the LC medium¹ in which the orientationally ordered materials are expected to exhibit surprising functionalities. Confinement is one of the many ways to induce phase behavior in LCs. Different types of geometric confinement, surface roughness, and particle–surface interaction can significantly affect the ordering properties of LC materials.^2,3 Although some routes that induce desirable phase transitions including applied external fields,^4–6 temperature changes,^7–9 coadsorption of mesogenic species¹⁰ and chemical modification of the surface¹¹ have been introduced, a new strategy is still necessary to achieve a better control of the spatial alignments of LCs in practical application.

Liquid crystal phases with a static boundary are well observed and characterized experimentally and theoretically, but the influence of a dynamic boundary on phase behavior is still less understood. We focus on the effects of a periodically active boundary on passive particles. With a periodic motion of the boundary, the system is dissipative, resembling fluid and granular material driven systems. Such systems driven out of the equilibrium states by changes in chemical or field-based stimuli have gained a lot of attention recently.^12–14 For example, the mechanism of how the dissipative self-assembly of a synthetic gelator uses chemical fuel as an energy source was studied by Boekhoven et al.¹⁵ Grzybowski et al. reported that dynamic patterns of millimeter-sized magnetic disks are driven by a rotating permanent magnet.¹⁶ Granular media under horizontal or vertical vibrations have also been extensively studied by experiment, theory, and computer simulations.^17–19 By controlling energy influx and energy dissipation, such driven systems can produce ordered structures, such as cluster phenomena, ordered separation, collective swirling motions,²⁰ giant number fluctuations,^17,21 phase transition,^22–24 and so on. These self-assembled structures resulting from non-equilibrium systems promise to be good candidates for novel functional materials.

Computer simulation studies based on proper modeling can provide guidance for experimental research and also offer a more direct means to study molecular systems, allowing the spatial behavior to be characterized in detail. In our paper, we investigate the structural transition confined in a simple periodically moving slit pore formed by two identical walls that are parallel to each other. Boundary vibration effectively changes the anchoring conditions instead of changing the effective interaction potential energy. It functions as the continuous input of energy in the monolayer by collisions between the boundary and ellipsoids. The energy is transmitted and dissipated by particles’ Brownian motion. A rich variety of configurations and the state phase diagram in the plane of amplitude and frequency are presented. We characterize these steady states by calculating local field parameters and energies. By varying the amplitude and frequency of the boundary, our study can provide a route to capture a range of ordered structures, which are expected to have potential applications in novel responsive materials.

Molecular simulation has been a useful method in studying the properties and dynamics of confined LCs. Here, we study a system of ellipsoid particles in an oscillating slit-like pore and probe the effect of the dynamic boundary on the formation of different structures and phase transitions by Langevin dynamics simulations. The Gay–Berne (GB) potential is a useful intermolecular potential and has been widely used for modeling of LC systems because it contains anisotropic repulsive and attractive interactions.^25–27 LCs are modeled by ellipsoids which interact via a generalized Gay–Berne potential. We focus on the simplest case in a two-dimensional system where the wall is a straight line. Each wall is composed of 51 spherical beads with a diameter σ₀, and their interactions are neglected. Systems containing N = 735 ellipsoid particles with an aspect ratio k = 3 are confined in a dynamic slab geometry with an initial wall separation L₀ = 50σ₀, with the walls located at y₀ = 0 and y_max = 50σ₀. The distance L between the two walls is changed periodically as L = L₀ − A + A cos 2π ft due to the boundary vibration, where A is the oscillating amplitude, and f is the frequency. The length of the simulation box in the x-direction is L_x = 50σ₀. The Cartesian coordinate system is set to be at the bottom wall, where x and y axes are along the horizontal and vertical directions, respectively. In addition, the vibration direction is along the y axis. We focus on relatively high bulk densities corresponding to area fractions of 70%, where the bulk phase is nematic. LC particles interact through the Gay–Berne potential,

where $u → ^ i$ and $u → ^ j$ are the unit vectors along the principal axis of LC molecules i and j, respectively. Here, $r → i j = r → i − r → j$ is the intermolecular vector distance, and $r → ^ i j = r → i j r i j$ is the unit vector along the vector $r → i j$⁠. The orientation-dependent range parameter $σ ( r → ^ i j , u → ^ i , u → ^ j )$ is given by

and the potential energy well-depth anisotropy is expressed as

with

and

where ϵ₀ and σ₀ are the energy and length units. The molecular anisotropy parameter χ depends on the aspect ratio $κ= σ e σ s$⁠, and σ_e and σ_s are the length of the major and minor axes of ellipsoids, respectively. χ′ is related to the ratio of the potential energy well depths, that is, κ′ = ϵ_s/ϵ_e, where ϵ_s and ϵ_e are the side-to-side and end-to-end configurations, respectively. The expressions of χ and χ′ are given by

Here, the parameters μ and υ are empirically determined dimensionless exponents and used to tune the shape of the anisotropic interaction. Each GB model is defined by a set of four parameters κ, κ′, μ, and ν, and the original GB parameters (3, 5, 2, and 1) are proposed by Gay and Berne.²⁸ 

Since the wall particles are spherical, their interactions with GB particles are analogous to GB–GB interactions by setting $u → ^ i =0$ in Eqs. (1)–(3).²⁹ Therefore, the range parameter and the potential energy well-depth anisotropy between the boundary beads and the ellipsoid particles become

and

We consider a nonequilibrium system of LCs confined in a dynamic slab pore and perform Langevin dynamics simulations using the LAMMPS package.³⁰ In our work, the motion of particles is described by the Langevin equations of motion. The friction coefficients for both translation and rotation are fixed at 10 to ensure slow diffusion of particles.³¹ All physical quantities are expressed in reduced units: σ₀ = 1.0, ϵ₀ = 1.0, and m₀ = 1.0; σ₀, ϵ₀, and m₀ are the units of length, energy, and mass, respectively. The reduced number density ρ is given in units of $1/ σ 0 2$⁠. Newton’s equations of motion for the force and torque acting on an ellipsoid particle are integrated using the Velocity–Verlet algorithm with a time step Δt = 0.0005τ, where $τ= ( σ 0 2 m 0 / ϵ 0 ) 1 2$⁠. It requires more time to reach the steady equilibration states, and a smaller time step is adopted for large amplitude or high frequency. The cutoff for the intermolecular interactions is chosen as r_c = 4.0σ₀. Periodic boundary conditions are imposed on the x axis. The trajectory and snapshots of particles are visualized by the OVITO tool. At the beginning, the ellipsoids are randomly distributed in the simulation box. All runs are evolved for a sufficiently long time so that the ellipsoids reach steady states in all cases. The collected data are averaged over 5000 steady configurations.

In order to characterize the structural properties of the particles, three local fields—the density ρ(y), the uniaxial order parameter S(y), and the tilt angle ψ(y)—are introduced.³² To obtain the three local quantities along the y-direction, the simulation box is divided into 100 equidistant bins along the y axis, and the local fields are calculated in each bin. The bin density is obtained by $ρ ( y ) = ⟨ N e ( y ) ⟩ L x δ y$⁠, where N_e(y) denotes the number of ellipsoids with their center of mass located in a bin. S(y) is defined as the positive eigenvalue of the local order matrix $Q α β ( y ) = N e ( y ) − 1 ⟨ ∑ i = 1 N [ 2 u α ( i ) u β ( i ) − δ α β ] ⟩$⁠, where u_α(i) is the αth component of the ith particle orientation vector $u → ( i ) = ( cos θ i , sin θ i )$⁠. S measures the degree of the orientational order of the phase. If S = 0, the system is in an isotropic state; S increases as the number of molecular axes aligning with the director increases; if S = 1, the system is perfectly uniaxially oriented. The corresponding eigenvector is the local director field $n →$⁠. The tilt angle ψ(y) denotes the angle between the local director and the x axis.

Figure 1 shows typical snapshots for a system with a static and a moving slit boundary. It is found that there are great differences between static and dynamic confinement. For the equilibrium case in static slit confinement displayed in Fig. 1(a), the particles are observed to be orientationally ordered and spatially disordered at this area fraction. The wall–particle interaction favors the particles that are parallel to the wall, and the system adopts a planar configuration. For dynamic confinement, the structure of the particles can be manipulated by the amplitude or frequency of the boundary. From Figs. 1(b)–1(d), it can be seen that the ellipsoids experience an N_par–S_per–S_cs phase transition at a fixed frequency f = 2.0 as the amplitude increases. At low frequency and small amplitude, the ellipsoids adopt a morphology, N_par, similar to the equilibrium nematic phase observed for molecules in static confinement. As the oscillating strength enhances, the orientation of the ellipsoids is changed from parallel to perpendicular to the wall. When the amplitude is higher, the core–shell structure S_cs appears where the core consists of densely packed and ordered particles surrounded by a shell made up of few particles in complete disorder. The appearance of the steady core–shell structure is a genetic phenomenon on phase separation in active matter, such as quorum-sensing run-and-tumble bacteria³³ and binary particles with diffusivity difference.³⁴ These dynamic states are essentially unchanged as long as the driving is maintained, so they are named steady structures. On increment of frequency with a fixed amplitude A = 4.0, the system of particles also undergoes different structural arrangements, as shown in Figs. 1(d)–1(f). The whole system tilts to the wall with f = 1 and A = 4.0, as shown in Fig. 1(e). The density of particles is very high, and the motion is greatly hindered. A multi-domain structure appears due to the anisotropic geometry of the ellipsoid. Although it is not obvious, disclination lines appear at the interfaces between the domains. The energy input increases due to acceleration of the boundary oscillating velocity, which makes the particles tend to disorder. First, the ellipsoids near the wall become disordered, and the core still possesses orientation and translation order [see Fig. 1(d)]; then, the whole system becomes disordered, characterized by random positions and orientations of ellipsoids with a larger amplitude and higher frequency, as shown in Fig. 1(f).

In order to characterize and distinguish the four structures, we consider three local fields, namely, density ρ, orientational order parameter S, and tilt angle ψ, along the y-direction, as shown in Fig. 2. At small amplitude, the density [Fig. 2(a) blue line] is constant except in a small region near the walls for the nematic structure N_par. The order parameter is approaching 0.9 [Fig. 2(b) blue line], and the particles slightly align, tilted by their long axes parallel to the walls, as displayed by the tilt angle [Fig. 2(c) blue line]. By increasing the amplitude, the particles aggregate inwardly in order to avoid collision with the boundary. In the S_per phase, ρ, S(y), and ψ [see Figs. 2(a)–2(c) black line] exhibit strong oscillatory behavior due to perpendicular layering of particles. The ellipsoids organize themselves into a smectic with layered structure, and well defined peaks arise compared to the other phases. The value of the local orientational order and the tilt angle is around 1.0 and $π 2$⁠, respectively, indicating a clear smectic order and a nearly perpendicular director to the walls. It can be noted that the distance between adjacent peaks is nearly 3σ₀ about the long axis of an ellipsoid molecule. The period of these oscillations also corresponds to the interlayer spacing of the ellipsoids. When the amplitude is further increased and S_cs forms, due to the molecules close to the walls becoming isotropic and the inner particles aligning themselves parallel to the boundary, the local order of particles decreases from the central region to the wall [Fig. 2(b) green line]. The density is uniform and high [Fig. 2(a) green line], and the tilt angle [Fig. 2(c) green line] is small in the core region. For the S_tilt structure formed at low frequency and large amplitude, the particles tilt more heavily toward the boundary, especially at the center. The density, orientational order, and tilt angle also oscillate [see Figs. 2(a)–2(c) red line], but the peaks are more irregular than those in the S_per phase.

To systematically investigate the effect of boundary motion on the steady structure, the phase diagram is depicted as shown in Fig. 3 in the plane of frequency f and amplitude A of the boundary motion. In the case of the dynamic boundary, the input energy due to the boundary motion plays an important role in determining the structure. Six distinct phases in the parameter range are identified, which are the near-equilibrium nematic phase, N_par, and smectic phase, S_par, parallel to the wall, the smectic phase perpendicular to the wall, S_per, the smectic phase tilted to the wall, S_tilt, the planarly ordered core particle and shell disordered S_cs structure, and the orientationally and translationally disordered D phase. When the boundary starts to oscillate with a small amplitude and low frequency, the near-equilibrium planar nematic and smectic phases appear. Further increasing A or f may lead to the appearance of out-of-equilibrium S_per, S_tilt, and S_cs steady structures. With a larger amplitude or higher frequency, the disordered state may occur because the boundary oscillation is so intense that the collective motion of particles is broken.

The tilt angle $ψ ¯$⁠, representing the angle between the director of the whole system and the x axis, is used to quantify the structural transition. If the $ψ ¯$ is 0, it means the whole system aligns parallel to the wall. When $ψ ¯$ is $π 2$⁠, it means the particles are perpendicular to the wall. Figures 4(a) and 4(b) provide the tilt angle $ψ ¯$⁠, $E ¯ k i n$⁠, and $E ¯ p o t$ by averaging 1000 periods of boundary oscillation to reflect the structural transition after the steady structure is achieved. It shows that the structure undergoes N_par → S_per → S_cs transition with the increase in the amplitude at a fixed frequency f = 2.0, as shown in Fig. 4(a), where $ψ ¯$ reveals abrupt jumps. The abrupt change in $ψ ¯$ suggests ellipsoid orientations ranging from planar to perpendicular and then to inner planar in reference to the boundary. $ψ ¯$ is slightly higher for S_cs than that for N_par due to a part of the particles near the boundary being orientationally isotropic. As the frequency is increased at a fixed amplitude A = 4.0, as shown in Fig. 4(b), the structure changes from S_per → S_tilt → S_cs → D. An abrupt change still can be found in the transition S_per → S_tilt and S_cs → D, while the transition S_tilt → S_cs is not apparent. Average kinetic and potential energies are also given in Figs. 4(a) and 4(b), which display a continuous increase. The increase in kinetic energy is obviously larger than the increase in potential energy, indicating that kinetic energy is dominant in structural transition.

In order to illustrate the major factor maintaining the steady structure from the energy point of view, we present the kinetic energy, potential energy, and total energy for the four steady states in Figs. 5(a)–5(d). The kinetic energy is $E k i n ′ ( t ) = ∑ i = 1 N ( 1 2 m v i 2 ( t ) + 1 2 I ω i 2 ( t ) )$⁠, where i denotes the ith particle, the potential energy $E p o t ′$(t) originates in the GB potential of the system, and the total energy is the sum of kinetic energy and potential energy. Periodical oscillations are found in all profiles, where work is carried out on the system as the boundary contracts and the system works externally as the boundary expands. The frequency of energy oscillation is consistent with the frequency of boundary motion, reflecting the characteristics of the system driven by the boundary oscillation. The percentage of potential energy is increased, as shown in Figs. 5(a)–5(d), because of the inward strengthened aggregation of the particles leading to close packing of ellipsoids. It is observed that the kinetic energy accounts for a large proportion of the total energy in all four profiles, which suggests kinetic energy has a major effect in maintaining the steady structure, and most of the work performed on the system is converted to kinetic energy.

Translational and rotational kinetic energy distributions per particle for four steady structures N_par, S_per, S_cs, and S_tilt are given in Figs. 6(a) and 6(b), which are used to characterize particle motion. The translational and rotational kinetic energy are calculated from $E t r a = ⟨ 1 n ∑ i = 1 n 1 2 m v i 2 ⟩$ and $E r o t = ⟨ 1 n ∑ i = 1 n 1 2 I ω i 2 ⟩$⁠, respectively. Here, n represents the number of particles in y → y + Δy, and ⟨⋯⟩ represents the ensemble average. The translational motion is the main way in particle collective motion compared to the rotational motion, as shown in Fig. 6. The active boundary works as an energy input in the monolayer by collisions between boundary and particles, which results in the increase in particle mobility. The ellipsoids have low translational speed in the whole region for N_par, except the several particles close to the boundary [Fig. 6(a) blue line]. As the amplitude increases, the particles in the central region still remain motionless, while ellipsoids outside gain increased translational kinetic energy for S_per [Fig. 6(a) black line]. In addition, the curve fluctuates due to the characteristic of the structure. The motion of inner particles is significantly increased, and the slow speed region decreases with a larger amplitude in the S_cs phase [see Fig. 6(a) green line]. The translational energy is nonlinearly increased from the center to the wall and then declines near the wall due to reduced particles. For the S_tilt structure formed at low frequency and large amplitude, the profile [see Fig. 6(a) red line] is close to a parabolic increment from the center to the boundary, and the particles inside are mainly immovable due to crowded ellipsoids with limited motion. The rotational kinetic energy curves [Fig. 6(b)] are similar to the translational energy curves for these structures. It can be concluded that translational speed is dominant in the steady structure and vibrating walls lead to inhomogeneous motion profiles along the y axis.

In this article, we have introduced a novel confining slit vibration system to study colloidal LCs by means of Langevin dynamics simulations. A tunable boundary is used so that the orientation of the particles can be varied while interaction between the particles and boundary remains constant. The moving boundary leads to different alignments compared with static bulk confinement. We have analyzed the structures by calculating the local fields for four steady structures N_par, S_per, S_cs, and S_tilt. These structural transitions are attributed to the particle aggregation and the increased mobility due to boundary oscillation. Furthermore, the oscillating boundary induces inhomogeneous mobility along the y axis. Our results reveal a strategy for new adaptable materials by designing dynamic organizing principles. We expect this non-equilibrium model to be helpful in addressing the ordered steady structure and the response to external perturbations of liquid crystalline systems applied in displays, sensors, and active materials.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This work was supported by the National Natural Science Foundation of China (Grant No. 51802051) and the Science and Technology Base and Talent Special Project of Guangxi Province (Grant No. 2018AD19098). We are grateful to the High Performance Computing Center (HPCC) of Nanjing University for performing the numerical calculations in this paper on its blade cluster system.