Structures, thermodynamics and dynamics of topological defects in Gay–Berne nematic liquid crystals

Abstract

Topological defects are a ubiquitous phenomenon across different physical systems. A better understanding of defects can be helpful in elucidating the physical behaviors of many real materials systems. In nematic liquid crystals, defects exhibit unique optical signatures and can segregate impurities, showing their promise as molecular carriers and nano-reactors. Continuum theory and simulations have been successfully applied to link static and dynamical behaviors of topological defects to the material constants of the underlying nematic. However, further evidence and molecular details are still lacking. Here we perform molecular dynamics simulations of Gay–Berne particles, a model nematic, to examine the molecular structures and dynamics of +1/2 defects in a thin-film nematic. Specifically, we measure the bend-to-splay ratio K₃/K₁ using two independent, indirect measurements, showing good agreement. Next, we study the annihilation event of a pair of ±1/2 defects, of which the trajectories are consistent with experiments and hydrodynamic simulations. We further examine the thermodynamics of defect annihilation in an NVE ensemble, leading us to correctly estimate the elastic modulus by using the energy conservation law. Finally, we explore effects of defect annihilation in regions of nonuniform temperature within these coarse-grained molecular models which cannot be analysed by existing continuum level simulations. We find that +1/2 defects tend to move toward hotter areas and their change of speed in a temperature gradient can be quantitatively understood through a term derived from the temperature dependence of the elastic modulus. As such, our work has provided molecular insights into structures and dynamics of topological defects, presented unique and accessible methods to measure elastic constants by inspecting defects, and proposed an alternative control parameter of defects using temperature gradient.

---

1 Introduction

Liquid crystals (LCs) represent a range of condensed matter phases that exhibit features of both simple liquids and crystalline solids.^1,2 The nematic phase, in which LC molecules display orientational ordering without positional ordering, is most studied thanks to its wide spectrum of applications, including display technology,³ sensors,^4,5 directed self-assembly,^6,7 and autonomous materials.⁸ The direction-dependent material properties of an LC are rooted in the anisotropic shape of its constituent molecules. Molecular dynamics (MD) simulations can therefore provide molecular-level insights into a variety of unique phenomena in LCs.^9–12 Examples include revealing molecular self-assembly structure on the interface of an LC nanodroplet,¹³ elucidating surface anchoring conditions for LC-liquid^14,15 and LC-solid boundaries,¹⁶ and resolving the molecular structure of topological defects.^17,18 Despite the above-mentioned advances of molecular models in understanding the equilibrium properties of LCs, dynamical behaviors of LCs are much less studied by molecular simulations.^9,19–21

A widely used coarse-grained model of low molecular weight LCs is the Gay–Berne (GB) model.^22–24 GB particles interact through a modified Lennard–Jones (LJ) potential accounting for particle shape anisotropy and interaction anisotropy.^22–24 The GB model can capture isotropic, nematic, smectic, and columnar phases in LCs.^25–31 Moreover, by carefully picking its parameters, the GB model can well simulate certain common LC molecules, e.g., 4-cyano-4′-pentylbiphenyl (5CB)³² and p-terphenyl.³³ The advantages of the GB model against other atomistic models are its affordable computational costs and its convenience in studying the genuine physics of LCs.³⁴ Previous works have been devoted to measuring material constants of GB particles by different methods. For example, orientational elastic moduli of GB particles have been measured using the direct correlation function method,³⁵ density of states,³⁶ free energy calculations,^34,37 and equilibrium orientational fluctuations.^38,39 Viscous coefficients of nematic GB systems were measured using equilibrium MD^40–46 and non-equilibrium MD simulations.^47–53 These fundamental studies have revealed microscopic structures and properties of LCs, and have paved the way toward more quantitative study of LC materials using molecular models.

Indeed, the GB model has been used to quantitatively elucidate a range of LC-related phenomena, such as anchoring effects,⁵⁴ wetting,⁵⁵ droplet shape,⁵⁶ nanofilm,^57–59 nanoconfinement effect,^60–72 nanoparticle effect,^73–77 boundary absorption of LC molecules,^16,78 and the phase behavior of GB–GB mixtures^79–83 and GB–LJ mixtures.^84,85 Other physics, such as the effect of external field^86,87 and the addition of charges^88,89 or dipoles^90–92 are also introduced in the GB model, showing interesting phenomena such as reorientation⁹¹ and the tilted phase.⁹² Recent GB particle simulation works have been extended to discotic,⁹³ biaxial,^94–101 deformable,¹⁰² or even active¹⁰³ LCs. These efforts have shed light on the structures of LCs at the microscopic scale, which is difficult to observe in experiments and inaccessible by continuum simulations.

In spite of the above research progress, the study of topological defects using the GB model is scarce. Topological defects are a consequence of broken symmetry in ordered systems. Because of its ubiquity, topological defects are important for understanding a wide variety of phenomena in various physical systems.¹⁰⁴ Topological defects in nematic LCs are regions where the orientational ordering of the molecules is frustrated. Point defects in two-dimensional (2D) nematics and wedge disclinations in three-dimensional nematics can be characterized by their winding number or topological charge, defined as k = α/2π, where α is the angle by which the director rotated (positive for counterclockwise rotation and negative for clockwise rotation) after a counterclockwise traversal of the Burgers circuit surrounding the defect.¹⁰⁵ Because of the nematic symmetry, the topological charge of a 2D point defect can be multiples of half integer. The self energy of such defect is proportional to the square of its topological charge.¹⁰⁵ Therefore, defects carrying the lowest topological charge ±1/2 are the stable ones in 2D nematics, which will be heavily studied here.

There is a recent interest in studying the structures and dynamics of topological defects in nematic LCs,^8,106–108 thanks to its emerging applications in defect-directed self-assembly,¹⁰⁹ material transport,¹¹⁰ and micro reactors. Interestingly, recent studies also showed that topological defects can reveal material properties of the LC. For instance, the morphology of a +1/2 defect is recently used to infer the elasticity of a lyotropic nematic LC;^111,112 it is believed that during the annihilation event of a pair of ±1/2 defects, the +1/2 defect moves faster than the −1/2 defect due to hydrodynamic effects.^111,113 Despite the success of continuum models that can reproduce the above-mentioned phenomena, molecular insights into these mechanisms are still lacking.

Here, we use MD simulations to investigate morphological and dynamical properties of topological defects in a nematic represented by GB particles. We first measure the shape of equilibrated +1/2 defects constructed by GB models with varying parameters and infer their bend-to-splay ratio K₃/K₁. Subsequently, we use the distance between two +1/2 defects in a quasi 2D nematic confined to a disk to infer K₃/K₁. The two independent, indirect measurements of the elasticity ratio agree well with previously published results. We next study the annihilation event of a pair of ±1/2 defects and find qualitative agreement with experiments and hydrodynamic simulations in terms of defect speed disparity and defect orientation dependence. We further extend our model to elucidate the thermodynamics of defects. Specifically, our micro-canonical ensemble (NVE) simulation gives a simple way of measuring the LC's elastic modulus using the energy conservation law. We also show that higher temperature leads to a longer time of defect annihilation. Finally, we demonstrate one special feature of MD simulation by considering the dynamics of a +1/2 defect under a spatial gradient of temperature. We find that its moving speed is accelerated along the positive gradient direction but decelerated in the opposite direction. This can be understood by considering the negative coupling between the temperature and the elastic modulus. Taken together, our MD simulations confirm that dynamical phenomena observed in macroscopic systems are also present at the nanoscale, which also agree with the continuum theory. Our simulations provide simple methods to infer elastic moduli of a nematic by inspecting defect structures and thermodynamics. Further, the prediction of temperature gradient effect on defect dynamics shows its promise in controlling defects for applications in, for example, defect-directed self-assembly, molecular transport, and nano-reactors.

2 Model and simulation details

2.1 Gay–Berne potential

The Gay–Berne (GB) potential is an anisotropic version of the LJ potential to mimick prolate and oblate spheroid molecules.^22–24 The interaction potential between two uniaxial GB particles can be written as follows:^22–24where û_i and û_j are unit vectors representing the orientations of particle i and j, respectively. r⃑_ij = r_ijr̂_ij is the center-to-center vector pointing from particle j to i. σ₀ is the diameter of the cross-section of the particle normal to its orientation (Fig. S1 in ESI,† Section S1). The ‘softness’ of the potential d_w allows for appropriate scaling of oblate spheroid particles, set to 0.345 for the κ = 0.345 oblate spheroid GB particles and 1 for all prolate spheroid particles. R approximates the surface-to-surface distance by

R = r_ij − σ(û_i,û_j,r̂_ij) + d_wσ₀.

The range parameter σ(û_i,û_j,r̂_ij) takes the following form
where the ratio χ = (κ² − 1)/(κ² + 1) characterizes the asphericity of the liquid crystal molecules, and κ represents the aspect ratio of the GB particle (Fig. S1 in ESI,† Section S1). The molecular anisotropy potential ε(û_i,û_j,r̂_ij) in eqn (1) is given by

ε(û_i,û_j,r̂_ij) = ε₀ε^ν₁(û_i,û_j)ε^μ₂(û_i,û_j,r̂_ij),

where ε₀ is the characteristic energy representing the energy well depth for cross configuration (Fig. S1 in ESI,† Section S1). ε₁ and ε₂ are defined as
and
Parameters ν and μ are dimensionless. Parameter , where κ′ = ε_ss/ε_ee is the ratio of the potential well depths corresponding to the side-to-side (ss) and end-to-end (ee) configuration, respectively (Fig. S1 in ESI,† Section S1). A cut-off distance r_cut is introduced to eqn (1) such that the potential U* = 0 when r_ij ≥ r_cut. A GB model can be fully determined by four parameters, denoted by GB(κ, κ′, μ, ν). In this paper, we use three kinds of GB particles, i.e., GB(4.4, 20, 1, 1), GB(3, 5, 1, 3) and GB(0.345, 0.2, 1, 2). Elementary units, i.e., energy scale ε₀, length scale σ₀, mass scale m₀, and the Boltzmann constant k_B, are set to 1. All units are multiples of these fundamental values, with distance x scaled as x* = x/σ₀, volume V as V* = V/σ₀³, energy U as U* = U/ε₀, temperature T as T* = k_BT/ε₀, density ρ as ρ* = Nσ₀³/V (N is the number of particles), pressure P as P* = Pσ₀³/ε₀, time t as t* = t(k_BT/m₀σ₀²)^1/2, elastic modulus K_i as = K_iσ₀/ε₀, and viscosity η as η* = ησ₀³/ε₀τ₀. In what follows, units are omitted if simulation units are used.

2.2 Simulation details

MD simulations were carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).¹¹⁴ Periodic boundary conditions are applied in all three directions, and the Nosé–Hoover thermostat is adopted to control the temperature in the isobaric-isothermal ensemble (NPT) and canonical ensemble (NVT). The time step is dt = 0.001 for GB(4.4, 20, 1, 1), GB(3, 5, 1, 3) and dt = 0.0001 for GB(0.345, 0.2, 1, 2) in simulation units for all investigations unless otherwise specified. The cut-off distance for GB(4.4, 20, 1, 1), GB(3, 5, 1, 3) and GB(0.345, 0.2, 1, 2) is r_cut = 6σ₀, 5σ₀ and 1.6σ₀, respectively. The NPT ensemble is set up using pressures extracted from NVT ensemble simulations, which is variable depending on the model and temperature used. The densities of GB(4.4, 20, 1, 1), GB(3, 5, 1, 3) and GB(0.345, 0.2, 1, 2) are ρ* = 0.1932, 0.3, and 2.36, respectively. To characterize nematic–isotropic phase transition, we use a cube as the simulation box (Fig. 1(a)). For the rest of the prolate GB systems, the z-dimension is much smaller than the other two dimensions and we treat them as quasi 2D systems. The simulation details of the parameter setting are summarized in Table 1.

Fig. 1 Morphology and structure of +1/2 defects in Gay–Berne nematics. (a) Phase behavior of a bulk nematic of GB(4.4, 20, 1, 1). (b) Scalar order parameter S as a function of temperature for three types of GB particles with comparison to the ref. 36. (c) GB(4.4, 20, 1, 1) particles confined to a disk with a homeotropic anchoring boundary condition. Defects morphology displayed by particle for GB(4.4, 20, 1, 1) (d), GB(3, 5, 1, 3) (g) and GB(0.345, 0.2, 1, 2) (i). r in (d) indicates the radial distance used to measure the morphology of +1/2 defects. Defects morphology displayed director field for GB(4.4, 20, 1, 1) (e), GB(3, 5, 1, 3) (h) and GB(0.345, 0.2, 1, 2) (k). The color bar indicates the order parameter S and dashed blue lines highlight the shape of the +1/2 defects. Two +1/2 defects observed in a disk region for GB(4.4, 20, 1, 1) (f), GB(3, 5, 1, 3) (i) and GB(0.345, 0.2, 1, 2) (l). Inset in (c): Colormap indicates the particle orientation for (a), (c), (d), (g) and (j).

Table 1 A summary of simulation parameters for different systems

GB(κ, κ′, μ, ν)	GB(4.4, 20, 1, 1)	GB(3, 5, 1, 3)	GB(0.345, 0.2, 1, 2)
(phase)	21.06 × 21.06 × 21.06	18.85 × 18.85 × 18.85	9.63 × 9.63 × 9.63
(disk)	(88–176) × 12	(60–120) × 8	(21.04–42.09) × 20
(annihilation)	143.44 × 143.44 × 28.17	163.55 × 163.55 × 22.36	63.44 × 63.44 × 15.16
(T gradient)	356.84 × 143.80 × 28.24	—	—
N (phase)	1805	2016	2106
N (disk)	15 540–59 220	7368–28 552	18 480–66 520
N (annihilation)	112 000	180 000	144 000
N (T gradient)	280 000	—	—

---

We use the Q-tensor formalism to calculate the scalar order parameter S and the director field in the nematic. The tensor representation is based on a spatial average of the dyadic product of individual molecular orientations. It is defined as^115,116

where N^r is the number of GB particles in the sampling region and I is the identity tensor. The local average orientation n̂ of the region is the eigenvector of Q associated with its largest eigenvalue S, which corresponds to the scalar order parameter in the local region quantifying its degree of nematic ordering. The size effect study of the sampling region shows that the scalar order parameter S is only marginally dependent on N^r (ESI,† Section S2). To measure the hydrodynamic velocity field in our GB system, we divide the simulation box into 15 × 15 regions in the xy plane, and use individual particle's time-dependent position vector to calculate a local region's average velocity v⃑. For the k-th region, the velocity vector v⃑_k is computed fromwhere i = 1, 2,…, denotes the i-th particle in the region, and x⃑_i(t*) represents the i-th particle's position vector at time t*. The time interval for the measurement is chosen to be Δt* = 200. Note that this hydrodynamic velocity vector v⃑ is different from the instantaneous velocity V⃑ of individual GB particles.

3 Results

We first study the nematic ordering of a bulk GB particle system at different temperatures using an NVT ensemble. For the temperature range we considered, nematic and isotropic phases are respectively detected below and above a critical temperature T_NI (Fig. 1(b)). Our measured scalar order parameter S as a function of temperature T* quantitatively agrees with the benchmark results reported in literature³⁶ for all three GB parameters, serving as a validation of our calculations (Fig. 1(b)).

Topological defects in a nematic LC are associated with high elastic energy region.¹⁰⁵ In continuum theory, the Frank–Oseen elastic energy density of a nematic LC is expressed in the following:¹¹⁷

where K₁, K₂, K₃, and K₂₄ are the splay, twist, bend, and saddle-splay modulus, respectively. In a 2D system, twist (K₂) and saddle-splay (K₂₄) terms are irrelevant. Therefore, the morphology of +1/2 defects is solely determined by the competition between splay (K₁) and bend (K₃) constants, which controls the in-plane distortions of the director field around the defect core.^111,112,118

To simulate the defect annihilation process in a GB system, we prepare a thin-film nematic containing a pair of ±1/2 defects separated by a distance d* = 16–70 and equilibrate these defects under an NPT ensemble (pressure obtained from NVT ensemble simulations) to reach the desired density for a short duration of t* = 5. Hereafter, we perform simulations in the NVT ensemble during which the two defects approach each other and eventually annihilate driven by the elastic force. We further measured the scalar order parameter S and the director n̂ around the core of the +1/2 defect based on eqn (2) well before the pair of defects are annihilated (Fig. 1(e), (h) and (k)), and the equilibrium S is different across the three types of GB particles. The low scalar order parameter areas (dark color) defined by S < 0.5 correspond to the defect cores and their sizes are (Fig. 1(e), (f), (h), (i), (k), (l) and ESI,† Section S3). The comet-like shape of the director field around the defect core is also slightly different among the three types of GB systems (Fig. 1(e), (h) and (k)). The defect associated with κ = 4.4 and 0.345 exhibits a V-like and a U-like shape, respectively. Whereas the defect associated with the intermediate aspect ratio κ = 3 shows an intermediate shape. To quantitatively characterize the shape of the +1/2 defect and to avoid boundary effect, we measure the azimuth dependence of the director orientation angle θ around the defect core at a fixed radial distance r* from the defect core before two defects are too close. This radial distance r* = 8, 8, 4 for κ = 4.4, 3, 0.345, respectively, which are chosen much shorter than the box dimension to minimize the influence of the other defect or the far field (Fig. 1(d), (g) and (j)). Note that the local Frank elastic constants close to defect core regions could deviate from their equilibrium values due to the lowering of the scalar order parameter S caused by the diverging elastic distortions. Therefore, r* is chosen to be considerably larger than .

To quantify the shape of +1/2 defect, we set a polar coordinate system represented by (r,ϕ) centered at the defect core and measure the orientation angle θ, defined as the angle between the y-axis and the director orientation n̂, as a function of the azimuth angle ϕ at the fixed radial distance r* (Fig. 2(a) inset). Our measurements indeed reveal the quantitative difference among the three GB systems we study here (Fig. 2(a)). By fitting our data with Frank elasticity theory,¹¹¹ we can further extract the elasticity ratio K₃/K₁. We find that the higher the aspect ratio κ is, the larger K₃/K₁ is. Moreover, κ > 1 and <1 lead to K₃/K₁ > 1 and <1, implying that the bend constant K₃ is larger (smaller) than the splay constant K₁ in a rod-(disk-)like system. These observations are consistent with our expectations; as past results have suggested that shape (rather than enthalpy or molecular identity) comprises most of the information necessary to capture liquid crystal ordering and elasticity.^36,119,120

Fig. 2 Two independent measurements of the elasticity ratio. (a) Morphology of +1/2 defects averaged from 40 data points for each GB particle type. Inset: Quantitative description of +1/2 defect morphology: θ is the angle between the y-axis and the director orientation; ϕ is the polar coordinate. (b) Defect distance measurement in disk simulations and theory, b is the extrapolation length. Error bars indicate standard deviation of 40 measurements; inset: d is the distance between two +1/2 defects in a disk with a homeotropic anchoring by fixing the GB particles (one layer near the boundary) along the radial direction.

We can also measure the elasticity ratio from an alternative system, namely a disk region consisting of two +1/2 defects. Specifically, we equilibrate a disk region of diameter D of GB particles with homeotropic anchoring by fixing the GB particles (one layer near the boundary) along the radial direction with time step dt = 0.0001. The separation distance between the two defects d is an outcome of the competition between splay and bend energy cost. Our MD simulations show that the average separation distance d̄ of the two defects is different among the three types of GB systems (Fig. 1(f), (i) and (l)). The higher κ is, the closer the two defects are. By minimizing the free energy of a nematic confined to a disk region using Ginzburg–Landau equation¹²¹ (ESI,† Section S4), we can theoretically find the equilibrium separation distance d of the two +1/2 defects with a known elastic constant ratio K₃/K₁ (Fig. 2(b)). The higher K₃/K₁ is, the more costly bend distortion is in the system, which will bring the two defects closer; otherwise, the two defects will move away from each other to suppress splay deformation. In the special case when K₃/K₁ = 1, our calculation shows d/D ∼ 0.667, consistent with what is reported in literature.^122,123

By pairing the elasticity ratio K₃/K₁ from the defect shape measurement and the defect separation distance d̄/D from the disk simulation, the different size data points (K₃/K₁,d̄/D) corresponding to the three different GB particle types collapse onto the theoretical curve (Fig. 2(b)). The slight deviation between data points and theoretical curve can be understood by considering effects of finite anchoring and unevenly distributed scalar order parameter (ESI,† Section S5). This remarkable result demonstrates that Frank elasticity theory still works in the nanoscale.

The elastic constant ratio reported from the density of state method³⁶ and our measurements are summarized in the following table:

Table 2 Thermodynamic quantities and elastic constants of the three GB systems

GB(κ, κ′, μ, ν)	GB(4.4, 20, 1, 1)	GB(3, 5, 1, 3)	GB(0.345, 0.2, 1, 2)
ρ*	0.1932	0.30	2.36
T*	3.00	2.80	2.60
T NI ^36	7.00	3.80	5.50
S ^36	0.74	0.82	0.73
K 1 ^36	2.43	5.00	6.00
K 3 ^36	6.07	8.32	4.05
K 3/K1^36	2.50	1.67	0.67
K 3/K1 (shape)	2.68	1.58	0.75
K 3/K1 (disk)	2.17–2.71	1.18–1.72	0.51–0.67

---

For both prolate and oblate (discotic) particles, three different measurements of K₃/K₁ differ by less than 12% (Table 2). The two methods provide a simple visual way to measure the elastic constant ratio in nematic LCs.

After studying the static properties of topological defects using the GB model, we turn to studying their dynamical behaviors. In what follows, we focus on GB(4.4, 20, 1, 1) system. A pair of ±1/2 defects in an otherwise uniform thin-film nematic is prepared to study their annihilation process. We consider two scenarios of defects configurations, namely a parallel scenario and a perpendicular scenario¹¹¹ (Fig. 3(a) and (b)). In a parallel (perpendicular) scenario, the auxiliary line connecting (dashed white line in Fig. 3(a) and (b)) the two defect cores is parallel (perpendicular) to the nematic far-field, the orientation of which is represented by the zoomed red and green particles in Fig. 3(a) and (b). In both scenarios, the two defects will approach each other to reduce the elastic energy of the system. The comet-like +1/2 defect will therefore move forward (backward) with respect to its head in the parallel (perpendicular) scenario (Fig. 3(a) and (b)). After the annihilation, the system returns to a uniform state.

Fig. 3 Two scenarios of defect annihilation events. (a) Sequential snapshots for the parallel scenario, the far field orientation (top left inset: the orientation of far field represented by the zoomed red particle) and the auxiliary line (dashed white line) connecting the two defects are parallel. Top right inset: Colormap indicates the particle orientation. (b) Sequential snapshots for the perpendicular scenario, the far field orientation (top left inset: the orientation of far field represented by the zoomed green particle) and the auxiliary line (dashed white line) connecting the two defects are perpendicular. (c) Defects position in x-axis as a function of time. The origin is set to the position of −1/2 defects in MD simulation; error bars indicate standard deviation of 10 independent simulation runs. (d) Defects position in x-axis as a function of time in LBM simulation with elastic constant ratio K₃/K₁ = 2.50. ξ and τ are the length and time unit in LBM simulations, respectively.

We have run 10 independent simulations starting from the same defect separation distance for each scenario, and measured the averaged defect trajectories shown in Fig. 3(c). For comparison, we have also performed hybrid lattice Boltzmann method (LBM) simulation using the same elastic constant ratio (K₃/K₁ = 2.50).¹²⁴ Our calculations show that the speed of the −1/2 defect is lower than that of the +1/2 defect (Fig. 3(c)), in line with the literature¹¹¹ and the hydrodynamic simulations (Fig. 3(d)), indicating that a hydrodynamic effect that causes defect speed disparity also exists in our nanoscale system.^111,113 We also notice that the two scenarios show different annihilation times. Specifically, the perpendicular scenario has faster annihilation dynamics than the parallel scenario (Fig. 3(c)), again consistent with the actin-based experiment¹¹¹ and the hydrodynamic simulations (Fig. 3(d)). Moreover, both MD and LBM simulations agree on the relative position of the annihilation point: the annihilation point of the two defects for the parallel scenario is closer to the original position of the −1/2 defect; whereas the annihilation point of the perpendicular scenario is slightly closer to the center point between the original positions of the two defects. Therefore, there is a semi-quantitative agreement between the two simulations in terms of the defect trajectories. No significant size effect on defect annihilation event is observed in our simulation (ESI,† Section S6). We further measure the velocity and the scalar order parameter field in the MD simulation (Fig. 4(a) and (c)) and contrast them to the hydrodynamic simulation results (Fig. 4(b) and (d)). The velocity vector field measured from the GB system based on eqn (3) is similar to that calculated from the hydrodynamic simulations (Fig. 4) and the literature.¹⁰⁸ During the annihilation, the prolate spheroid GB particles rotate and translate collectively to accommodate the re-arrangement of the director field, leading to a hydrodynamic flow (Fig. 4(a) and (c)). For both scenarios, two vortices of flow are seen distributed on the two sides of the auxiliary line connecting the two defects, which push the +1/2 defect to move towards the −1/2 defect (Supplementary Movie, ESI†). Whereas there is no significant large-scale flow near the −1/2 defect. This backflow effect explains the observation that the +1/2 defect moves faster than the −1/2 defect in GB systems (Fig. 3(c)). The annihilation time difference between the two scenarios can be understood by considering the viscous coefficients associated with the backflow. In the parallel scenario, the director field between the two defects is perpendicular to the flow. For the perpendicular case, however, the director field between the two defects is aligned with the flow. Therefore, the relevant Miesowicz viscosity¹²⁵ in parallel and perpendicular scenarios is η₂ and η₁, respectively. Because η₂ > η₁ in GB(4.4, 20, 1, 1),⁴⁶ the hydrodynamic flow is stronger in perpendicular scenario than in parallel scenario, therefore explaining the annihilation time difference between the two scenarios.

Fig. 4 Velocity field comparison between MD and LBM simulation for two defect annihilation scenarios. Snapshots for the parallel scenario (a) and perpendicular scenario (c) in MD simulation. Snapshots for parallel scenario (b) and perpendicular scenario (d) in LBM simulation with elastic constant ratio K₃/K₁ = 2.50. Blue arrows indicate the velocity vector. Color indicates scalar order parameter S. τ is the time unit in LBM simulations.

Although we averaged the velocity in each region, the velocity fluctuations in MD simulations are more significant than that in LBM simulations. The color in Fig. 4 of two methods look differently. That is because thermal fluctuations of the order parameter are absent in an LBM simulation, which is based on a continuum model. In MD simulations, thermal fluctuations are naturally included and are reflected by the non-uniform color in the images. For equilibrium simulations, this noise can be reduced by doing longer time average. For dynamic simulation in which defects move and annihilate, the instantaneous scalar order parameter will be noisy. In MD simulations, the displacement of the defects in the y-direction fluctuates, making the vortex center and vertex shapes change during the annihilation (Fig. 4(a) and (c)). Moreover, there is a diverging flow at the final stage of the annihilation process in the LBM simulation, which is not observed in the MD simulation. Nevertheless, an acceleration of the −1/2 defect speeds is found in MD simulation, consistent with the fact that the elastic force diverges as their separation distance approaches zero, driving a strong flow in the final stage of the annihilation process (Fig. 3(c)).

To estimate the Ericksen number of the GB model in the defect annihilation event, we have measured the rotational viscosity using rotating field method (ESI,† Section S7.1). The characteristic velocity of the two defects v̄ = 0.036σ₀/τ₀ (Fig. 3(c)), the core size of defects during annihilation ξ = 2, and choose the mean elastic constant K* = 4.25.³⁶ This gives rise to v̄ξ/K* = 0.224, comparing to Er = 0.219 in the LBM simulation. Therefore, the viscous effect is considerable but less important than the elastic force for this range of the Ericksen number. This is also consistent with the fact that the hydrodynamic effect is present in our GB systems. Note that, similar to the concept that one can estimate orientational elastic moduli by inspecting defect shapes, one can also estimate viscous coefficients by inspecting defect velocities. Indeed, this method yields , which is in the same order of magnitude of the rotating field measurement (ESI,† Section S7.2).

We next study the thermodynamics of defect annihilations using different ensembles. We first consider an NVE ensemble and set the initial temperature to T* = 3. During annihilation, we detect a steady increase of the temperature despite strong fluctuations in the measurements (Fig. 5(a)). After annihilation, temperature fluctuates with respect to a higher mean value with T* ≈ 3.005 (Fig. 5(a)). The increase of the thermal energy associated with the temperature rise allows us to estimate the elastic modulus of the nematic. Specifically, the elastic energy corresponding to the initial state of the two defects separated by a distance d* = 70 can be approximated by a elastic energy formula under one-constant approximation:¹⁰⁵

where k₁ = 1/2 and k₂ = −1/2 are the topological charge of the two defects, respectively, K* is the mean elastic modulus, and is the characteristic size of the defect core. After annihilation, this elastic energy is dissipated into heat in the NVE ensemble. The thermal energy change in the simulation is:

E_thermal = Nk_BΔT*.

We can use the relationship E_thermal = E_el to estimate the elastic constant,
The estimation of elastic constant is K* = 4.61, agreeing well with (K₁ + K₃)/2 = 4.25 reported in literature.³⁶ We also find that when the distance d* is comparable to the defect core size, the elastic constant measurement using the energy conservation law will become inaccurate (ESI,† Section S8). This provides a convenient method to estimate the order of magnitude of the elastic modulus of the nematic in MD simulations. This estimate, if combined with the visual measurement of the elasticity ratio, can be used to determine the absolute values of both K₁ and K₃. This idea of applying thermodynamic frameworks to study defect annihilation dynamics is a fertile area for future works.

Fig. 5 Defect annihilation simulations at different ensembles and temperatures. (a) Temperature change during defect annihilation in the NVE ensemble. (b) Pressure change during defect annihilation in the NVT ensemble. (c) Volume change during defect annihilation in the NPT ensemble. (d) Defect annihilation time as a function of defects annihilation point in x-axis in different ensembles (NVT, NVE, and NPT) for both parallel and perpendicular scenarios. (e) Defects position in x-axis as a function of time with different temperatures. Error bars indicate standard deviation of 10 independent simulation runs.

To ensure that our calculations are independent of the choice of the ensemble, we simulate defect annihilations in otherwise equivalent NPT, NVT and NVE ensembles using GB(4.4, 20, 1, 1). A change in thermodynamic variables during defect annihilation is also observed. In the NVT ensemble, the measured pressure P has dropped by ∼0.35%; whereas, in the NPT ensemble, the measured volume has shrunk by ∼0.17%. Note that the presence of defects can distort the nematic and lead to less efficient packing of the GB particles, resulting in higher pressure for the fixed volume system and larger volume for the fixed pressure system. In all the three ensembles we consider here, the fluctuations of thermodynamic variables, 〈Δx〉/〈x〉 (x = T, P, or V) is in the order of 10⁻³, comparing to . This shows that our nanoscale system is within the thermodynamic limit (ESI,† Section S9). To compare defect annihilation dynamics between different ensembles, we conduct NVT, NVE, and NPT ensemble simulations with the same initial temperature and density, and compare their annihilation dynamics in terms of annihilation time t* and the location of annihilation point x* (Fig. 5(d)). There is no significant difference between different ensembles within statistical uncertainty, confirming that the GB system is indeed in the thermodynamic limit and different ensembles are indifferent.

We further investigate the temperature dependence of defect annihilation processes using the NVT ensemble. The effect of temperature has two folds. On the one hand, viscous coefficients are usually lower at higher temperature,⁴⁶ which can facilitate hydrodynamic flows and help accelerate defect annihilation. On the other hand, elastic moduli are lower at higher temperature,³⁶ which will reduce the driving force of the defects, slowing down the annihilation process. Our MD simulations can be used to ascertain the relative importance of these two competing effects. Our measured defect trajectories at three different temperatures show that higher temperature leads to a longer time of annihilation (Fig. 5(e)), implying that elasticity is more important than the viscous effect in terms of annihilation dynamics in the GB system considered here. This is again consistent with the estimated Ericksen number Er = 0.224 < 1 of the system at which elastic effect wins over viscous effect.

Finally, we study the effect of temperature gradient on defect dynamics. It is advantageous to use MD simulations against continuum simulations to study non-uniform temperature systems. The temperature dependence of material constants such as elastic moduli and viscous coefficients has to be implemented ad hoc in continuum simulations. In molecular models, however, they emerge naturally. In what follows, we use the GB model to investigate how a +1/2 defect moves under a temperature gradient. To this end, we consider a pair of ±1/2 defects separated in the x-direction along which the local temperature varies linearly from at the center to at the two edges of the box, where the periodic boundary condition requires that their temperatures are equal (Fig. 6(a) and (b)). If , the +1/2 defect facing towards the center is under a positive temperature gradient (Fig. 6(a)); if otherwise, the +1/2 defect is under a negative temperature gradient (Fig. 6(b)). The temperature range of the system bounded by and are within the nematic phase. Our simulations show that the +1/2 defect under a positive temperature gradient spends less time to move to the center of the box than under a negative temperature gradient (Fig. 6(c) and (d)). We also notice that GB particles tend to rotate out of the plane near the center of the box, a similar behavior found in uniform temperature systems (Fig. 3). This phenomenon is more severe for the positive temperature gradient scenario in which the box center is hotter, where stronger thermal fluctuations may promote the out-of-plane rotations of the GB particles.

Fig. 6 Temperature gradient effect for defect motion. Diagram of temperature distribution along the x-axis of positive (a) and negative (b) temperature gradient. Spatiotemporal kimograph of defects motion for positive (c) and nagative (d) gradient scenario. The color bar in (b) indicates the temperature magnitude. Inset in (d): colormap indicates the particle orientation.

We further quantify the change of moving speed of +1/2 defects under different temperature gradients. We subtract the defect velocity v_x by v_x0 measured from a simulation of a uniform temperature and plot Δv_x = v_x − v_x0 against in Fig. 7(a). The results show that the more positive (negative) the temperature gradient is, the higher (lower) the defect speed will be.

Fig. 7 Temperature gradient effect on +1/2 defect motion. (a) Change of +1/2 defect speed Δv_x as a function of temperature gradient ∂T*/∂x*. (b) Temperature gradient (along y-direction) effects for defect motion. The color bar indicates the temperature magnitude. Inset: colormap indicates the particle orientation. Error bars indicate standard deviation of 40 measurements.

We can understand the above behavior using a continuum theory. The elastic energy of the system is under one-constant approximation,¹⁰⁵ where K* is a function of temperature T*.³⁶ The force acting on the +1/2 defect can be derived by differentiating the elastic energy , where the first term contributes to the change of defect speed Δv_x and the second term corresponds to v_x0 for a uniform temperature system. Therefore, the theory shows that there is a linear relation between Δv_x and ∂T*/∂x* and the proportionality is , comparing to 3.95 measured from a linear fit to the simulation data (Fig. 7(a)). This quantitative agreement again demonstrates that the molecular system can be understood from a continuum point of view. The above predictions can be further tested in future experiments. Note that a recent experiment of a colloid dispersed in nematic 5CB has reported a similar phenomenon that the colloid tends to move toward the hotter region.¹²⁶ This is due to a similar physical reason that the elastic modulus of the hotter regions is lower, which provides a force to push the colloid to the hot zone to reduce the elastic distortion costs incurred by the presence of the colloid, an opposite direction when dispersed in an isotropic medium under the same temperature gradient.¹²⁶ To further support our proposed mechanism, we perform a separate simulation in which a temperature gradient is along the transverse direction (y-direction), where the periodic boundary condition is turned off to suppress any global flow (Fig. 7(b)). The results also show that the +1/2 defect tends to move towards the hotter region, which provides an additional support of these theoretical arguments.

4 Discussion and conclusion

In this work, we have performed molecular dynamics simulations of topological defects in a thin-film Gay–Berne nematic. We first use two independent, indirect methods to measure the bend-to-splay ratio K₃/K₁ for three types of GB particles, including two different prolate particles and one discotic particle. In the first method, we inspect the shape of the +1/2 defect and fit it with a continuum theory. In the second method, we measure the distance of two +1/2 defects of a nematic confined to a disk region with homeotropic anchoring. By comparing to the continuum theory, K₃/K₁ measured by the two methods agree with each other and are within ∼10% difference with the literature.³⁶ This provides a direct molecular evidence of the Frank elasticity theory. We next study the spontaneous annihilation process of a pair of ±1/2 defects in an otherwise uniform nematic of GB(4.4, 20, 1, 1). The defect trajectories indicate that the +1/2 defect moves faster than the −1/2 defect, and the configuration of the director field can also dictate the annihilation speed. Both features are thought to originate from hydrodynamic effects of the nematic. This agreement with experiments and continuum simulations confirms the hydrodynamic theory of nematics. We can also extract the flow field of the GB system during defect annihilation and compare favorably with that calculated from hydrodynamic simulations. We further examine the thermodynamics of the system before and after defect annihilation. The NVE ensemble simulation allows us to correctly estimate the elastic modulus by equating the elastic energy stored in the defect-bearing state and the thermal energy acquired after defect annihilation. Combined with the visual methods of measuring the elasticity ratio, it is possible to use the information of defect structures and dissipated energy during defect annihilation to determine the absolute values of K₁ and K₃. By varying temperature, simulations show that defect annihilation is faster for lower temperature, underscoring the importance of the elastic effect, which becomes stronger at lower temperature. This is consistent with our estimated Ericksen number of the system Er = O(0.1) at which elastic force wins over viscous force.

Temperature gradients can lead to spatial gradients of elastic constants, viscous coefficients, and diffusion constants. Landau–de Gennes theory can be used in the quasi-static limit and when thermal noise is not considered. In general, a hydrodynamic model accounting for elastic constant gradients, back-flow effects, and thermal diffusion effects has to be developed to fully describe the thermophoresis phenomenon of defects. In addition to the complexity of this model, the temperature dependence of these material constants is unknown, with limited measurements available from which to estimate their magnitude, therefore requiring they be set ad hoc in the model as well. In the molecular model, however, the temperature dependency of these constants naturally arises and no ad hoc assumption is required, and all the physics are included. Therefore, the molecular model provides a convenient platform onto which to theoretically study the temperature gradient effect. Finally, we simulate defect motion in a region of nonuniform temperature. Our calculations show that +1/2 defects tend to move to hotter area. Using a simple theory based on the temperature dependence of the elastic modulus, we are able to quantitatively explain the change of defect speed due to local temperature gradient. There is a recent interest in studying the transport phenomena of topological defects in nematic liquid crystals. For example, certain types of defects can self propel in active nematics.¹²³ Because defects tend to segregate foreign molecules and particles in nematic liquid crystals,^6,109,127 controlled motion of defects motion can lead to applications in, for example, defect-carried material transport and micro reactors. In active nematics, internal stresses have to be maintained to mobilize these defects. Therefore, temperature gradient provides an alternative method to mobilize defects locally and in demand.

Put together, our work presents a compendium of new investigations into the physics of topological defects in nematic liquid crystals; we have systematically compared MD and continuum simulations in terms of defects structures, thermodynamics, and annihilation dynamics. The study shows that there is a quantitative agreement between the two physical models of distinct length scales, including the measured elastic constant ratio, annihilation dynamics, and temperature gradient effect on defect motion. Moreover, the measurement of viscoelastic constants of a material is often-times cumbersome and inevitably involves the application of external fields. In this work, we propose passive methods to measure these constants, namely through inspecting defect structure, defect separation distance, or measuring temperature change during an annihilation event. The temperature gradient effect provides an alternative way to control defect dynamics using heat, and presents a novel prediction which should stimulate future experimental interest, particularly as thermophoretic effects are likely to be present in the burgeoning field of active liquid crystals.

Author contributions

Y. H. performed MD simulations and analyzed the data. W. W. performed LBM simulations and analyzed the data. R. Z. designed, supervised the study, and performed continuum calculations. Y. H. and R. Z. interpreted the data and wrote the first draft. All authors contributed to the revision of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge support from the Research Grants Council of Hong Kong SAR through grant no. 16300221 and the HKUST Central High Performance Computing Cluster (HPC3). JKW acknowledges support from the US National Science Foundation, Award No. DMR-1751988. The authors also thank Tiezheng Qian, Xinyu Wang and Shengzhu Yi for fruitful discussions. The authors also thank the Referees for their inspiring and helpful suggestions.

Notes and references

1. P.-G. De Gennes and J. Prost, The physics of liquid crystals, Oxford University Press, 1993.
2. P. Oswald and P. Pieranski, Nematic and cholesteric liquid crystals: concepts and physical properties illustrated by experiments, CRC Press, 2005.
3. R. H. Chen, Liquid crystal displays: fundamental physics and technology, John Wiley & Sons, 2011.
4. I.-H. Lin, D. S. Miller, P. J. Bertics, C. J. Murphy, J. J. de Pablo and N. L. Abbott, Science, 2011, 332, 1297–1300.
5. M. Sadati, A. I. Apik, J. C. Armas-Perez, J. Martinez-Gonzalez, J. P. Hernandez-Ortiz, N. L. Abbott and J. J. de Pablo, Adv. Funct. Mater., 2015, 25, 6050–6060.
6. J. K. Whitmer, X. Wang, F. Mondiot, D. S. Miller, N. L. Abbott and J. J. de Pablo, Phys. Rev. Lett., 2013, 111, 1–5.
7. M. Rahimi, T. F. Roberts, J. C. Armas-Pérez, X. Wang, E. Bukusoglu, N. L. Abbott and J. J. de Pablo, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 5297–5302.
8. R. Zhang, A. Mozaffari and J. J. de Pablo, Nat. Rev. Mater., 2021, 6, 437–453.
9. C. Care and D. Cleaver, Rep. Prog. Phys., 2005, 68, 2665.
10. R. Berardi, L. Muccioli, S. Orlandi, M. Ricci and C. Zannoni, J. Phys.: Condens. Matter, 2008, 20, 463101.
11. C. Zannoni, Liq. Cryst., 2018, 45, 1880–1893.
12. M. P. Allen, Mol. Phys., 2019, 117, 2391–2417.
13. J. A. Moreno-Razo, E. J. Sambriski, N. L. Abbott, J. P. Hernández-Ortiz and J. J. de Pablo, Nature, 2012, 485, 86–89.
14. M. Sadati, H. Ramezani-Dakhel, W. Bu, E. Sevgen, Z. Liang, C. Erol, M. Rahimi, N. Taheri Qazvini, B. Lin and N. L. Abbott, et al. , J. Am. Chem. Soc., 2017, 139, 3841–3850.
15. H. Ramezani-Dakhel, M. Sadati, M. Rahimi, A. Ramirez-Hernandez, B. Roux and J. J. de Pablo, J. Chem. Theory Comput., 2017, 13, 237–244.
16. V. Palermo, F. Biscarini and C. Zannoni, Phys. Rev. E, 1998, 57, 2519–2522.
17. M. Rahimi, H. Ramezani-Dakhel, R. Zhang, A. Ramirez-Hernandez, N. L. Abbott and J. J. de Pablo, Nat. Commun., 2017, 8, 1–8.
18. O. D. Lavrentovich, P. Pasini, C. Zannoni and S. Zumer, Defects in liquid crystals: Computer simulations, theory and experiments, Springer Science & Business Media, 2001, vol. 43.
19. P. Pasini and C. Zannoni, Advances in the computer simulations of liquid crystals, Springer Science & Business Media, 1999, vol. 545.
20. D. Cheung, S. Clark and M. R. Wilson, Chem. Phys. Lett., 2002, 356, 140–146.
21. Z. Ran, P. Zeng-Hui, L. Yong-Gang, Z. Zhi-Gang and X. Li, Chin. Phys. B, 2009, 18, 4380.
22. B. J. Berne and P. Pechukas, J. Chem. Phys., 1972, 56, 4213–4216.
23. J. Gay and B. Berne, J. Chem. Phys., 1981, 74, 3316–3319.
24. P. Prybytak, W. Frith and D. Cleaver, Interface Focus, 2012, 2, 651–657.
25. E. de Miguel, L. F. Rull, M. K. Chalam and K. E. Gubbins, Mol. Phys., 1991, 74, 405–424.
26. E. de Miguel, E. M. Del Rio, J. T. Brown and M. P. Allen, J. Chem. Phys., 1996, 105, 4234–4249.
27. M. A. Bates and G. R. Luckhurst, J. Chem. Phys., 1999, 110, 7087–7108.
28. V. V. Ginzburg, M. A. Glaser and N. A. Clark, Liq. Cryst., 1997, 23, 227–234.
29. E. de Miguel and C. Vega, J. Chem. Phys., 2002, 117, 6313–6322.
30. D. Caprion, L. Bellier-Castella and J.-P. Ryckaert, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 041703.
31. J. A. Moreno-Razo, O. Cienegas-Caceres, E. Díaz-Herrera and J. Quintana, AIP Conf. Proc., 2008, 979, 120–129.
32. H. Fukunaga, J.-I. Takimoto, T. Aoyagi, T. Shoji, F. Sawa and M. Doi, Mol. Cryst. Liq. Cryst. Sci. Technol., Sect. A, 2001, 365, 739–746.
33. G. Luckhurst and P. Simmonds, Mol. Phys., 1993, 80, 233–252.
34. A. A. Joshi, J. K. Whitmer, O. Guzmán, N. L. Abbott and J. J. de Pablo, Soft Matter, 2014, 10, 882–893.
35. A. Avazpour, Liq. Cryst., 2012, 39, 1491–1497.
36. H. Sidky and J. K. Whitmer, Liq. Cryst., 2016, 43, 2285–2299.
37. J. Stelzer, L. Longa and H.-R. Trebin, J. Chem. Phys., 1995, 103, 3098–3107.
38. M. P. Allen, M. A. Warren, M. R. Wilson, A. Sauron and W. Smith, J. Chem. Phys., 1996, 105, 2850–2858.
39. A. Humpert and M. P. Allen, Mol. Phys., 2015, 113, 2680–2692.
40. A. Smondyrev, G. B. Loriot and R. A. Pelcovits, Phys. Rev. Lett., 1995, 75, 2340.
41. S. Sarman, J. Chem. Phys., 1998, 108, 7909–7916.
42. A. Cuetos, J. M. Ilnytskyi and M. R. Wilson, Mol. Phys., 2002, 100, 3839–3845.
43. S. Sarman and A. Laaksonen, Chem. Phys. Lett., 2009, 479, 47–51.
44. S. Sarman and A. Laaksonen, J. Chem. Phys., 2009, 131, 144904.
45. G. R. Luckhurst and K. Satoh, J. Chem. Phys., 2010, 132, 184903.
46. K. Satoh, Mol. Cryst. Liq. Cryst., 2015, 615, 78–88.
47. N. Mori, J. Morimoto and K. Nakamura, Mol. Cryst. Liq. Cryst. Sci. Technol., Sect. A, 1997, 307, 1–15.
48. L. Bennett and S. Hess, Proceedings of the 26th Freiburger Arbeitstagung Flussigkristalle, 1997.
49. L. Bennett, S. Hess, C. Pereira Borgmeyer and T. Weider, Int. J. Thermophys., 1998, 19, 1143–1153.
50. L. Bennett and S. Hess, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 5561.
51. N. Mori, R. Semura and K. Nakamura, Mol. Cryst. Liq. Cryst. Sci. Technol., Sect. A, 2001, 367, 445–453.
52. C. Wu, T. Qian and P. Zhang, Liq. Cryst., 2007, 34, 1175–1184.
53. S. Sarman, Y.-L. Wang and A. Laaksonen, Phys. Chem. Chem. Phys., 2015, 17, 16615–16623.
54. J. A. Moreno-Razo, E. J. Sambriski, G. M. Koenig, E. Díaz-Herrera, N. L. Abbott and J. J. de Pablo, Soft Matter, 2011, 7, 6828–6835.
55. D. Vanzo, M. Ricci, R. Berardi and C. Zannoni, Soft Matter, 2016, 12, 1610–1620.
56. D. Vanzo, M. Ricci, R. Berardi and C. Zannoni, Soft Matter, 2012, 8, 11790–11800.
57. A. P. Emerson, S. Faetti and C. Zannoni, Chem. Phys. Lett., 1997, 271, 241–246.
58. S. J. Mills, C. M. Care, M. P. Neal and D. J. Cleaver, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 58, 3284.
59. S. J. Mills, C. M. Care, M. P. Neal, M. R. Wilson, M. P. Allen and D. J. Cleaver, J. Mol. Liq., 2000, 85, 185–195.
60. G. D. Wall and D. J. Cleaver, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 56, 4306.
61. R. E. Webster, N. J. Mottram and D. J. Cleaver, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 021706.
62. R. Latham and D. J. Cleaver, Chem. Phys. Lett., 2000, 330, 7–14.
63. E. Cañeda-Guzmán, J. A. Moreno-Razo, E. Díaz-Herrera and E. J. Sambriski, Mol. Phys., 2014, 112, 1149–1159.
64. D. Salgado-Blanco, C. I. Mendoza, M. A. Chávez-Rojo, J. A. Moreno-Razo and E. Díaz-Herrera, Soft Matter, 2018, 14, 2846–2859.
65. D. Salgado-Blanco, E. Daz-Herrera and C. I. Mendoza, J. Phys.: Condens. Matter, 2019, 31, 105101.
66. D. Caprion, Eur. Phys. J. E: Soft Matter Biol. Phys., 2009, 28, 305–313.
67. Q. Ji, R. Lefort and D. Morineau, Chem. Phys. Lett., 2009, 478, 161–165.
68. Q. Ji, R. Lefort, R. Busselez and D. Morineau, J. Chem. Phys., 2009, 130, 234501.
69. Q. Ji, R. Lefort, A. Ghoufi and D. Morineau, Chem. Phys. Lett., 2009, 482, 234–238.
70. R. Busselez, C. V. Cerclier, M. Ndao, A. Ghoufi, R. Lefort and D. Morineau, J. Chem. Phys., 2014, 141, 134902.
71. J. Karjalainen, J. Lintuvuori, V.-V. Telkki, P. Lantto and J. Vaara, Phys. Chem. Chem. Phys., 2013, 15, 14047–14057.
72. S. I. Hernández, J. A. Moreno-Razo, A. Ramírez-Hernández, E. Díaz-Herrera, J. P. Hernández-Ortiz and J. J. de Pablo, Soft Matter, 2012, 8, 1443–1450.
73. B. T. Gettelfinger, J. A. Moreno-Razo, G. M. Koenig, J. P. Hernández-Ortiz, N. L. Abbott and J. J. de Pablo, Soft Matter, 2010, 6, 896–901.
74. T. Stieger, M. Schoen and M. G. Mazza, J. Chem. Phys., 2014, 140, 1–10.
75. T. Stieger, S. Püschel-Schlotthauer, M. Schoen and M. G. Mazza, Mol. Phys., 2016, 114, 259–275.
76. O. Guzmán, E. B. Kim, S. Grollau, N. L. Abbott and J. J. de Pablo, Phys. Rev. Lett., 2003, 91, 1–4.
77. A. D. González-Martínez, M. A. Chávez-Rojo, E. J. Sambriski and J. A. Moreno-Razo, RSC Adv., 2019, 9, 33413–33427.
78. G. D. Wall and D. J. Cleaver, Mol. Phys., 2003, 101, 1105–1112.
79. O. Cienega-Cacerez, C. García-Alcántara, J. A. Moreno-Razo, E. Díaz-Herrera and E. J. Sambriski, Soft Matter, 2016, 12, 1295–1312.
80. R. A. Bemrose and C. M. Care, Mol. Phys., 1997, 90, 625–636.
81. S. J. Mills, C. M. Care, M. P. Neal and D. J. Cleaver, Mol. Cryst. Liq. Cryst. Sci. Technol., Sect. A, 1999, 330, 1667–1674.
82. J. A. Moreno-Razo, E. Díaz-Herrera and S. H. Klapp, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 19–21.
83. R. Berardi and C. Zannoni, Soft Matter, 2012, 8, 2017–2025.
84. D. Antypov and D. J. Cleaver, J. Phys.: Condens. Matter, 2004, 16, S1887.
85. R. Berardi, A. Costantini, L. Muccioli, S. Orlandi and C. Zannoni, J. Chem. Phys., 2007, 126, 044905.
86. R. Berardi, S. Orlandi and C. Zannoni, Phys. Chem. Chem. Phys., 2000, 2, 2933–2942.
87. E. E. Burnell, R. Berardi, R. T. Syvitski and C. Zannoni, Chem. Phys. Lett., 2000, 331, 455–464.
88. S. Orlandi, L. Muccioli, M. Ricci, R. Berardi and C. Zannoni, Chem. Cent. J., 2007, 1, 1–13.
89. M. Lamarra, L. Muccioli, S. Orlandi and C. Zannoni, Phys. Chem. Chem. Phys., 2012, 14, 5368–5375.
90. R. Berardi, S. Orlandi and C. Zannoni, J. Chem. Soc., Faraday Trans., 1997, 93, 1493–1496.
91. R. Berardi, S. Orlandi and C. Zannoni, Int. J. Mod. Phys. C, 1999, 10, 477–484.
92. R. Berardi, S. Orlandi and C. Zannoni, Mol. Cryst. Liq. Cryst., 2003, 394, 141–151.
93. O. Cienega-Cacerez, J. A. Moreno-Razo, E. Díaz-Herrera and E. J. Sambriski, Soft Matter, 2014, 10, 3171–3182.
94. R. Berardi, C. Fava and C. Zannoni, Chem. Phys. Lett., 1995, 236, 462–468.
95. D. J. Cleaver, C. M. Care, M. P. Allen and M. P. Neal, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 54, 559.
96. R. Berardi, C. Fava and C. Zannoni, Chem. Phys. Lett., 1998, 297, 8–14.
97. R. Berardi and C. Zannoni, J. Chem. Phys., 2000, 113, 5971–5979.
98. R. Berardi, M. Cecchini and C. Zannoni, J. Chem. Phys., 2003, 119, 9933–9946.
99. R. Berardi, J. S. Lintuvuori, M. R. Wilson and C. Zannoni, J. Chem. Phys., 2011, 135, 134119.
100. S. Orlandi, L. Muccioli and R. Berardi, Liq. Cryst., 2018, 45, 2400–2415.
101. L. Querciagrossa, M. Ricci, R. Berardi and C. Zannoni, Phys. Chem. Chem. Phys., 2017, 19, 2383–2391.
102. L. Muccioli and C. Zannoni, Chem. Phys. Lett., 2006, 423, 1–6.
103. R. F. de Souza, S. Zaccheroni, M. Ricci and C. Zannoni, J. Mol. Liq., 2022, 352, 118692.
104. J. J. Sandford O’Neill, P. S. Salter, M. J. Booth, S. J. Elston and S. M. Morris, Nat. Commun., 2020, 11, 1–8.
105. M. Kleman and O. D. Lavrentovich, Soft matter physics: an introduction, Springer, 2003.
106. J. Brugués, J. Ignés-Mullol, J. Casademunt and F. Sagués, Phys. Rev. Lett., 2008, 100, 037801.
107. J. Galanis, R. Nossal, W. Losert and D. Harries, Phys. Rev. Lett., 2010, 105, 168001.
108. X. Tang and J. V. Selinger, Soft Matter, 2019, 15, 587–601.
109. X. Wang, D. S. Miller, E. Bukusoglu, J. J. De Pablo and N. L. Abbott, Nat. Mater., 2016, 15, 106–112.
110. R. Zhang, S. A. Redford, P. V. Ruijgrok, N. Kumar, A. Mozaffari, S. Zemsky, A. R. Dinner, V. Vitelli, Z. Bryant and M. L. Gardel, et al. , Nat. Mater., 2021, 20, 875–882.
111. R. Zhang, N. Kumar, J. L. Ross, M. L. Gardel and J. J. de Pablo, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, E124–E133.
112. S. Zhou, S. V. Shiyanovskii, H.-S. Park and O. D. Lavrentovich, Nat. Commun., 2017, 8, 1–7.
113. G. Tóth, C. Denniston and J. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 051705.
114. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
115. K. Schiele and S. Trimper, Phys. Status Solidi B, 1983, 118, 267–274.
116. I. W. Stewart, The static and dynamic continuum theory of liquid crystals: a mathematical introduction, CRC Press, 2019.
117. F. C. Frank, Discuss. Faraday Soc., 1958, 25, 19–28.
118. S. D. Hudson and E. L. Thomas, Phys. Rev. Lett., 1989, 62, 1993.
119. H. Sidky, J. J. de Pablo and J. K. Whitmer, Phys. Rev. Lett., 2018, 120, 107801.
120. J. Shi, H. Sidky and J. K. Whitmer, Mol. Syst. Des. Eng., 2020, 5, 1131–1136.
121. M. Ravnik and S. Žumer, Liq. Cryst., 2009, 36, 1201–1214.
122. V. Vitelli and D. R. Nelson, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 051105.
123. A. Mozaffari, R. Zhang, N. Atzin and J. de Pablo, Phys. Rev. Lett., 2021, 126, 227801.
124. R. Zhang, T. Roberts, I. S. Aranson and J. J. de Pablo, J. Chem. Phys., 2016, 144, 084905.
125. M. Miesowicz, Nature, 1946, 158, 27.
126. J. Kołacz, A. Konya, R. L. Selinger and Q.-H. Wei, Soft Matter, 2020, 16, 1989–1995.
127. X. Wang, Y.-K. Kim, E. Bukusoglu, B. Zhang, D. S. Miller and N. L. Abbott, Phys. Rev. Lett., 2016, 116, 147801.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2sm01178f

---