Molecular simulations of entangled defect structures around nanoparticles in nematic liquid crystals

ABSTRACT

We investigate the defect structures forming around two nanoparticles in a Gay–Berne nematic liquid crystal using molecular simulations. For small separations, disclinations entangle both particles forming the figure of eight, the figure of omega and the figure of theta. These defect structures are similar in shape and occur with a comparable frequency to micron-sized particles studied in experiments. The simulations reveal fast transitions from one defect structure to another suggesting that particles of nanometre size cannot be bound together effectively. We identify the ‘three-ring’ structure observed in previous molecular simulations as a superposition of the different entangled and non-entangled states over time and conclude that it is not itself a stable defect structure.

Graphical Abstract

KEYWORDS:

• Molecular-simulation
• defects
• nematic
• disclination
• algorithmic classification

1. Introduction

Spherical particle inclusions in liquid crystals distort the orientational order of the surrounding molecules inducing topological defects. Defect regions are defined by the absence of orientational order and a significant biaxiality [1]. Such defects arise due to the competition between liquid crystal molecules trying to align along the average molecular direction, called the director, and trying to fulfil the surface anchoring condition of the nanoparticle. For small particles or particles confined to a thin nematic cell, the defect region forms a defect ring near the equator of the particle with respect to the director, commonly referred to as a Saturn ring [2]. Particles surrounded by a Saturn ring show quadrupolar interactions due to the symmetry of the director field surrounding them; hence, such particles are called quadrupoles. For two quadrupoles in close vicinity, entangled defects were found to spontaneously arise when the surrounding nematic is distorted. Using laser tweezers, which allow easy manipulation with great precision, a range of reproducible entangled objects has been found [3–9]. Here, ‘entangled’ means that both particles are surrounded by a single defect line, commonly referred to as a disclination line. The three main structures comprise the figure of eight (), the figure of omega () and the figure of theta (). Schematics are shown in Figure 1. Micrographs of these entangled defect structures observed in experiments can be found in Figure 4 of Tkalec and Muševič [10]. For the figure of eight defect, a single disclination line winds around the particles in the manner shown in Figure 1(a). The figure of omega consists of one disclination line that surrounds both particles near the equator, entering into the space between the particles to describe a shape similar to the Greek letter . The figure of theta consists of a disclination line surrounding both particles and an additional defect loop in the plane in between the particles. Note the spontaneous symmetry breaking for the first two of these structures due to their chirality. The two states are degenerate and hence left-handed and right-handed structures occur with the same frequency.

Figure 1. (Colour online) Schematic of two nanoparticles (grey) in nematic host (not shown) entangled by a disclination line (dark red) forming the (a) figure of eight , (b) the figure of omega and (c) the figure of theta .

Figure 2. (Colour online) Schematic of simulation box. (a) Typical snapshot of Gay–Berne molecules near a Lennard-Jones wall applied near the z boundaries indicated in grey. Gay–Berne molecules were colour coded according to their orientation with respect to the director. (b) Two-dimensional view of the simulation box (grey dashed region) in the x–z plane. The black box indicates the area used to simulate an infinitely long one-dimensional chain of nanoparticles (blue). The box length along y was unchanged. Molecular snapshot shows a slice of the new simulation box with the nanoparticles at its centre.

Figure 3. (Colour online) Disclination lines around two nanoparticles (grey). Blue small dots correspond to data points with below the threshold and red large points correspond to points calculated using the ‘sphere-tracking’ (see description in text).

Figure 4. (Colour online) Schematic of two nanoparticles (grey) and the surrounding disclination lines (black) in the x–y plane. Red regions show the starting areas for the flood fill to distinguish different defect structures. Depending on which red regions are connected the defect type was determined.

These entangled defects were observed in experiments and predicted numerically for micron-sized particles using phenomenological mean-field calculations based on Landau–de Gennes free energy minimisation [7,11]. At first glance, these entangled structures appear to disobey the restriction that the net topological charge within the liquid crystal must be zero; however, Čopar and Žumer [8] resolved this apparent contradiction with the introduction of a self-linking number.

A closer look at the different defect structures reveals that entangled defects differ only within a tetrahedral region. Theoretically, by rotating the tetrahedron, one entangled state can be transformed into another [8]. In practice, this is equivalent to cutting and reconnecting the disclination lines, whereas in experiments, this rewiring is achieved by locally applying laser tweezers to heat the region around the tetrahedron.

Molecular simulations of two quadrupoles in close vicinity have shown a three-ring structure for particle inclusions of nanometre size [3,9,12]. This structure is similar to the figure of ; however, the inner loop connects with the outer loop at two nodes. Two possible explanations have been proposed as to why the structure differs from the ones observed in experiments. Earlier studies declared it a transient state [4]. It was later proposed [13] that instead the size of the particles is crucial and that the three-ring structure appears for particles of nanometre scale, whereas micron-sized particles form only the three different entangled defects (, and ). In this paper, we propose an alternative explanation for the three-ring structure.

The system of two entangled particles can be thought of as a building block for more complex systems. Recent work has shown one-dimensional [5] as well as two-dimensional [14–16] aggregates, where the disclination line winds around all the particles. Entanglement has also been observed for microspheres and microfibres [17]. When bringing many particles into close vicinity, knotted structures can be created [16,18,19]. By applying the laser tweezers in unknotted regions, additional knots can be created and vice versa. Wood et al. [20] showed that liquid crystals with a high number of nanoparticle inclusions can be used to synthesise a soft solid, which shows high rigidity due to a network of entangled defect lines. These materials have important features for potential use in biosensors [21].

Entangled and knotted systems can be manipulated not only by the use of laser tweezers but also by using strong local electric fields, hydrodynamic flow, temperature changes and the use of different boundary conditions and confined geometries. This allows the creation and modification of defects in liquid crystals in a rather controllable way and hence a variety of complex systems can be designed.

In this paper, we outline a model to efficiently simulate entangled nanoparticles. Simulation details are given in Section 2. In Section 3, we present how the different defect types are automatically categorised. The results, in particular the time scales and dynamics of transitions from one entangled defect structure to another, are discussed in Section 4. In addition, the origin of the three-ring structure, the only structure observed in previous molecular simulations, is identified. Conclusions are drawn in Section 5.

2. Model and simulation details

The entanglement was studied by simulating several systems with nanoparticle inclusions of various sizes and separations using the molecular dynamics package LAMMPS [22]. The nematic host was simulated using the well-known soft Gay–Berne potential (see Appendix A), a coarse-grained single-site potential that represents the interaction energies between two elongated particles. We chose a length-to-width ratio of . The remaining Gay–Berne coefficients were set to , and and a potential cut-off of was applied. Note that reduced units were used throughout by setting and the molecular mass leading to a basic unit of time . The moment of inertia was set to assuming uniform mass distribution. The initial director orientation was chosen to lie along the z direction.

Spherical particle inclusions of radius were placed in the box, interacting with the liquid crystal molecules via a purely repulsive Lennard-Jones potential (see Appendix B). Their positions were fixed throughout the simulation. The system was simulated for particles separated along x with a surface-to-surface separation of . Five different systems of nematics with nanoparticle inclusions of radii , and were studied. In all cases, box dimensions were chosen to give a bulk density , and the temperature was chosen to be , which, for this system, gives a nematic liquid crystal phase. The presence of the nanoparticles is expected to affect the local pressure tensor only in the immediate vicinity of the particle surfaces, not in the bulk nematic phase, and therefore, we assume that the usual isotropic box scaling algorithm is sufficient. Simulation details are given in the following subsections.

2.1. Nanoparticles of radius

Two spherical particle inclusions of radius were placed in the box, separated along x, with as well as . A cubic simulation box of length with periodic boundaries was used. The large simulation box is necessary to avoid any interference between nanoparticles through the boundaries. The system was equilibrated over 1.6 × 10⁶ time steps using an ensemble, followed by a production run, also in the ensemble, of 2.5 × 10⁵ time steps with a time step . To simulate the state point described above, the corresponding pressure was set to . Molecular positions and orientations were stored every 500 time steps. The total number of Gay–Berne molecules was 512,000.

2.2. Nanoparticles of radius

The simulation was repeated with particles of size separated along x with as well as . Most other simulation details were unchanged. Two modifications were made to reduce computational cost.

The first modification consists of the introduction of flat walls at the z boundaries. A simple Lennard-Jones 12-6 function of the z coordinate, relative to the wall, was utilised with and a potential cut-off of . This has the advantage that can be chosen to be smaller than in the system described above, without the danger of long-range interactions between images of nanoparticles across the boundaries; this type of boundary also stabilises the bulk director in the z direction. must be sufficiently large to accommodate the density and order parameter fluctuations near the walls, ensuring bulk behaviour in the centre of the simulation box. A system snapshot of Gay–Berne molecules near a Lennard-Jones wall can be seen in Figure 2(a). One can clearly see layering near the wall that disappears quickly further away from the wall.

The second modification consists of the reduction of the simulation box length to (so for and for ) and the placement of one nanoparticle at (), illustrated in Figure 2(b). With periodic boundary conditions along x, this represents an infinite chain of nanoparticles along the x direction ( remains large enough to avoid periodic self-interaction). This allows the reduction of the total number of Gay–Berne molecules by over 50%; choosing in both cases, with for and for , gave a bulk density . The system was equilibrated over 7.5 × 10⁵ time steps () using an ensemble, followed by a production run (NVT ensemble) of 4 × 10⁵ time steps. Molecular positions and orientations were stored every 500 time steps.

2.3. Nanoparticles of radius

In this simulation, the nanoparticle radius was increased to . Other system parameters were chosen as for the system with . No wall potential was applied and instead periodic boundaries were used across z. The total number of Gay–Berne molecules was with a moment of inertia of . The box dimensions were chosen to be and ; hence, the separation between neighbouring nanoparticles was . The system was equilibrated over 8 × 10⁵ time steps () using an NVT ensemble, followed by a production run (NVE ensemble) of 6.5 × 10⁵ time steps. Molecular positions and orientations were stored every 500 time steps.

3. Automated characterisation of defect types

To identify defect regions for the system snapshots stored, a weighted order tensor was calculated following the approach suggested by Callan-Jones et al. [23]. This quasi-continuous tensor is given by

(1)

where is the number of molecules in the spherical sampling volume with radius centred at , is the position vector of particle i and is the component of its orientation vector with . is a weighting function with the constraint

(2)

We choose a cubic b-spline, which is a piecewise continuous cubic polynomial approximation to a Gaussian function [24]:

(3)

where . is zero if is greater than the radius . was chosen to be , which corresponds to having roughly 30 Gay–Berne molecules inside the sampling volume . A scaling was imposed afterwards on w such that Equation (2) was obeyed. The eigenvalues of D are labelled . The uniaxial order of a grid point is defined as its alignment to the local nematic director field and is given by the Westin metric , which is zero in the isotropic region and unity in perfectly ordered regions. Similarly, a Westin metric can be defined, which is directly proportional to the biaxiality of the local director field. Within the bulk, and undergo small variations, whereas significant changes indicate defect regions. We identified as an appropriate threshold to determine defect core regions. By choosing a suitable threshold for , we were able to identify almost identical defect regions. For convenience, we choose to focus on from here onwards.

For each system snapshot, the weighted order tensor D was calculated on a regular 3D grid with a spacing of and defect regions were evaluated, which allowed the subsequent determination of the defect structure. Due to the high number of snapshots, this process was automated (https://github.com/sausages/sausages) and the defects classified as follows.

All grid points of the defect core, corresponding to , were identified and stored in a set ; all points with values above the threshold were discarded. We did not further distinguish between different values of .

To separate the points of into disclination lines, a queue-based flood-fill algorithm was used, a pseudocode representation of which is given in Algorithm 1. To begin, an arbitrary starting point is removed from and added to an empty queue. The algorithm then consists of moving the first point p from the queue to a list and moving any neighbours of p from to the end of the queue, in any order. This step is repeated until the queue is empty, at which point, the list corresponds to a discrete contiguous volume of points, which we identify as a disclination line. A new arbitrary starting point is then moved from to the newly empty queue to form the seed of a new list , and the process continues until is empty and so all points have been assigned to a disclination line.

Algorithm 1 Flood-fill based algorithm to separate points into disclination lines

Set i to 1

while is not empty do

Choose a point at random

Remove from

Add to queue Q

while Q is not empty do

Remove the first point from Q

Add to set

for each neighbour of do

Remove from

Append to Q

end for

Increment i

end while

Note that since only defect regions were stored, a point may not have neighbours in all six directions. To efficiently implement the flood-fill algorithm for such an irregular array, a linked-list connectivity graph was constructed such that each point stores references to its neighbours. Since disclination lines can only exist in the form of closed loops (or lines across periodic boundaries), each defect line was checked for this criterion. For coarser grids, noise artefacts in the data sometimes caused hairline breaks in the disclination lines, causing them to have ‘endpoints’. These endpoints could be identified by examining the neighbour-to-neighbour distance between points, estimated using the all-pairs shortest-path Floyd–Warshall algorithm [25] to find the two points separated by the longest path. Two endpoints in close proximity, and with no other endpoint nearby, could be automatically joined and their lists merged.

While the Floyd–Warshall algorithm gives a deterministic measure of distance along the disclination line, it scales with the cube of the number of points P and would be too slow for larger systems. The sparsity of the point-connectivity graph would make a per-point Dijkstra approach more promising in these cases, with scaling [25]. Coarse graining could also be used to reduce the effective P.

The length of a disclination line is of particular interest because it is proportional to its associated energy. Most of our disclination lines were well behaved, allowing the use of a ballistic ‘sphere-tracking’ approach which is much faster than connectivity-based approaches. An arbitrary starting point is chosen, and a second is chosen nearby. We define the unit vector connecting these two points as our (arbitrary) initial direction of travel and move a short distance along it. The average position of the defect line near the destination was calculated by averaging the positions of the 10 nearest points, which was enough to encapsulate the cross section of the line for our grid size. The new direction of travel is set to be the vector connecting the previous position and this new average, and this direction is followed from the most recent average to arrive at the new destination. The process continues until the initial starting point is reached. The length was estimated by summing the distances between neighbouring averaged positions. This method proved to be accurate and reliable. In Figure 3, the red points indicated the average positions calculated using this method for an example snapshot.

To accurately distinguish between the different defect types, the connections between the red regions in Figure 4 were evaluated using the flood-fill algorithm. Here, the boundaries were treated as fixed. The connections in conjunction with the number of disclinations and their respective lengths allow us to determine the defect structure. For the figure of eight, we further distinguished the direction of the twist, depending on which path crosses over the other. In over 99.9% of the snapshots, this automated analysis was conclusive; inconclusive exceptions were inspected visually.

4. Results and discussion

In Figure 5, the different defect structures over time are shown for all five systems. In total, five different structures were observed: separated Saturn rings, figure of eight, figure of omega, figure of theta as well as an intermediate structure, where the two distant parts of the disclination are linked. With the exception of the intermediate defect structure, all these defect structures were also observed in experiments and using Landau–de Gennes minimisation for micron-sized particles [26]. For each of these structures, a typical example is shown in Figure 6.

Figure 5. Different defect structures shown over time for all five simulations. Particle radius and surface-to-surface separation are given by R and , respectively. The defect types are labelled as following: intermediate defect structure (linked), figure of omega (), figure of theta (), two separate Saturn rings (2 SR) and figure of eight (). Subscripts indicate right- and left-handed twist.

Figure 6. (Colour online) Defect structures (red) observed for two nanoparticles (grey) of radius in close proximity (). Some are shown from two angles for clarity. (a) Two separate Saturn rings (2 SR), (b) figure of eight (), (c) figure of omega (), (d) figure of theta () and (e) intermediate defect structure (Linked).

The only non-entangled structure found was two nanoparticles surrounded by a Saturn ring each shown in Figure 6(a). The Saturn rings are strongly bent away from each other to minimise the distortion of the director field, which minimises free energy. The bending effect was also observed for Landau–de Gennes free energy minimisations [26]. Figure 6(b–d) shows typical snapshots of the figure of eight, the figure of omega and the figure of theta from two different angles. The shape of the defect line is in good agreement with the observations for micron-sized particles in experiments. Figure 6(e) shows an intermediate structure, where two distant segments of the disclination line are linked. For all five structures, thermal fluctuations of the disclinations are visible, since time averaging was avoided.

For nanoparticles with radius , the defect structure transitions frequently between different structures and no structure persists for longer than . When comparing the results for the two different separations, one can see that transitions are more frequent for the smaller separation, suggesting a smaller barrier to interconversion. It appears, at least for the time simulated, that two separate Saturn rings are the most frequent structure. The figure of eight can be seen occasionally, whereas the figures of theta and omega are very rare. The exact frequencies are given in Table 1 for all five simulations. The intermediate linked structure was observed frequently. This is interesting as this structure is not observed in experiments, suggesting that it is unstable, which in turn suggests that the particle size does indeed have a significant impact on the energy barriers of the entanglement.

Table 1. Frequency of observations of different entangled defect structures for a different particle radii and separations given in %.

	Defect structure frequency (in %)
	2 SR					Linked
2.94	31.7	5.9	16.8	1.0	2.0	42.6
10.30	71.9	3.6	13.4	0.0	0.0	11.2
5.0	24.3	22.2	9.0	0.0	5.2	39.2
10.0	99.88	0.0	0.0	0.0	0.0	0.12
10.0	94.5	3.6	0.0	0.0	0.0	1.8

and are the particle radii and their surface-to-surface separation, respectively. Notation for different defect structures as introduced in Figure 6.

For particles with radius and a surface-to-surface separation of , entangled defect structures formed, and the frequency of transitions is very similar to that for the smaller particles. With the separation increased to , no transitions were observed throughout the production run and the only structure seen was the two separate Saturn rings. For particles with radius and a surface-to-surface separation of , two Saturn rings were observed most frequently. For a short period, , the figure of eight with a right-hand twist formed.

To summarise, it appears that two separate Saturn rings are the most stable structure for the nanoparticle sizes studied here. When the particles are sufficiently close, entangled and intermediate structures were observed with frequent transitions between them. The figure of eight was observed to be more stable than the figure of omega and the figure of theta was the least stable.

At this point, it has to be clarified that results presented in Table 1 are indicators only. Much longer simulation run times would be necessary to give accurate numbers. For instance, the linked defect structure depends on the resolution grid size chosen when calculating the weighted order tensor, i.e. a structure may appear to be linked for low resolution but could be classified otherwise if a higher resolution was used. We chose a resolution of expecting this to be smaller than the size of the defect core and hence giving reasonably accurate results. Furthermore, Table 1 shows that the frequencies of the figure of eight defect with left-handed and right-handed twist are dissimilar. This underlines that the simulation run times used were insufficient to obtain accurate results. Further investigation, for examples of the free energies associated with the different entangled structures, would have to take this into account; however, it is beyond the scope of this paper.

These concerns aside, our simulations clearly indicate that two separate Saturn rings are more stable than entangled structures and that the figure of eight is more stable than the figure of omega, which in turn is more stable than the figure of theta. This is in quantitative agreement with experimental results for micron-sized particles in strong confinement [5]: two separate quadrupoles are the only stable structure, observed in 48% of all laser tweezers manipulations. For entangled structures, the figure of eight was found most frequently (36%), followed by the figure of omega (13%). The figure of theta was observed rarely (3%). However, the binding energies seem to be much smaller for nanoparticles, indicated by the frequent transitions. By contrast, micron-sized particles’ binding energies were calculated to be an order of magnitude stronger than the unbound pair and several thousand times stronger than particles dispersed in water [5].

We tried to analyse the relation between the length of the disclination line, the occurrence of a transition and the state the defect structure is in. However, we found no significant relation. This is partially due to the high fluctuations in length, where the line stretched or shrank by a few within 500 time steps and also due to the inaccuracy of the ‘sphere-tracking’ measurement that only provides estimates within . Nonetheless, it would be an interesting feature to study for more stable entangled defects.

Finally, we address the three-ring structure that was observed in previous molecular simulations [3,9,12]. In Figure 7, the isosurface corresponding to an order parameter is plotted for the time-averaged results for the entire production run for and . The time-averaged disclination line forms the three-ring structure with two nodes previously observed. It appears that it is a product of the frequent transitions and not a stable defect structure itself. We do not observe such a structure in any instantaneous snapshot. It appears that the visualisation technique here is important: the common approach using small bins is inadequate to capture the fast fluctuations. Using a weighted order tensor including several neighbouring molecules on a fine grid has proven itself as a very accurate method to locate defect regions.

Figure 7. (Colour online) Disclination line (red) corresponding to time-averaged over entire production run for and .

Our results are in good agreement with calculations by Araki et al. [4,15], who observed the figure of eight for nanoparticles using numerical methods based on nematodynamics. In addition, we have observed the figure of omega and figure of theta for the first time in molecular simulations.

The results suggest that the entanglement for nanoparticles of sizes studied here is too weak to observe stable entangled defects, given that the defects frequently transition between different states and often form intermediate structures which we expect to be highly unstable. Future work should address larger particle inclusions to generate more stable defect structures; this could be achieved using the techniques described in Section 2. Whether a molecular simulation could access the longer length scales necessary to observe transitions between stable entangled defects remains to be seen.

5. Conclusions

Molecular simulations were successfully used to simulate defects around two spherical nanoparticles in close vicinity in a nematic host. In our exploratory simulations, we studied three different radii and several different surface-to-surface separations. Five different defect structures formed: well-separated Saturn rings, the figure of eight, the figure of omega, the figure of theta and an intermediate structure. To our knowledge, this is the first observation of the figure of omega and figure of theta for nanoparticles. The observation of intermediate structures was due to the small size of the particle inclusions as well as simulations allowing access to much smaller time scales than experiments. All defect structures observed were qualitatively similar to observations made for micron-sized particles; however, for the particle sizes studied here, the transitions were very fast and none of the entangled structures persisted for more than a few hundred time steps. This suggests that very small particles cannot be effectively bound together by entangled lines and that thermal energies are higher than the energy barriers between different entangled defect structures. Our analysis reveals that the three-ring structure observed in previous molecular simulations is solely an artefact of time-averaged superpositions of the different entangled states and not a stable defect structure itself.

To further explore the phenomenon of entangled nanoparticles in molecular simulations, we suggest a significant increase in particle size. We expect that the transitions will occur less frequently for larger particles and hence the dynamics of the system could be studied. However, one has to keep in mind that this would require much longer simulation times to study transitions. At this point, we have to stress that our simulations were pushing the current limits of computing resources available. Therefore, it is likely that rare event simulation methods (see e.g. Refs. [27,28]) will have to be applied in larger simulations to explore all types of entanglement.

Acknowledgements

The computer facilities provided by Warwick University Centre for Scientific Computing and by the ARCHER UK National Supercomputing Service are gratefully acknowledged as well as the support from the Engineering and Physical Sciences Research Council.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The computer facilities provided by Warwick University Centre for Scientific Computing and by the ARCHER UK National Supercomputing Service are gratefully acknowledged as well as the support from the Engineering and Physical Sciences Research Council.

Appendix A.

The Gay–Berne potential

The potential originally suggested by Gay and Berne [29] is widely used to simulate liquid crystals. It can be regarded as a shifted, anisotropic Lennard-Jones potential, i.e. it depends on the relative orientation of the particles as well as their separation. For identical uniaxial particles, it can be written as [30,31]:

(A1)

where

(A2)

and are unit vectors along the principal axes of the two particles i and j, while is the vector connecting their centres of mass and . is a parameter representing the width of the particle. is the orientation-dependent range parameter

where

Here, is given by , where is the length-to-width ratio of the particle. The strength anisotropy function used in Equation (A1) is given by

(A3)

is the well-depth parameter determining the overall strength of the potential, while and are two adjustable exponents which allow considerable flexibility in defining a family of Gay–Berne potentials. and are given by

where

Here,

(A4)

where is the ratio of well depths for the side-to-side configuration, , and the end-to-end configuration, , of two molecules.

Appendix B.

Homeotropic surface anchoring

Instead of using a specific anchoring potential, a simple variation of the standard Lennard-Jones (LJ) 12-6 potential is used. Here, the anchoring is entirely induced by the packing effects of the liquid crystal particles near the surface of the nanoparticle. For the homeotropic surface anchoring, the Gay–Berne (GB) molecules are allowed to penetrate the surface of the nanoparticle. To prevent GB molecules from entirely entering the particle, a shifted purely repulsive LJ interaction potential

(B1)

is used. is an energy parameter chosen to be unity and is given by

(B2)

Here, is the vector connecting the positions of the nanoparticle and the GB molecule and is the corresponding unit vector. is a size parameter and defined as the smallest diameter of the GB molecule; in this system, . is the distance of closest approach between the GB molecule and the nanoparticle, set to

(B3)

where is the radius of the spherical nanoparticle. For this interaction potential, the potential cut-off is chosen to be .