Defect-mediated colloidal interactions in a nematic-phase discotic solvent

Abstract

Interactions between colloidal inclusions dispersed in a nematic discotic liquid-crystalline solvent were investigated for different solute–solvent coupling conditions. The solvent was treated at the level of Gay–Berne discogens. Colloidal inclusions were coupled to the solvent with a generalized sphere-ellipsoid interaction potential. Energy strengths were varied to promote either homeotropic or planar mesogenic anchoring. Colloid–colloid interactions were modeled using a soft, excluded-volume contribution. Single-colloid and colloid-pair samples were evolved with Molecular Dynamics simulations. Equilibrium trajectories were used to characterize structural and dynamical properties of topological defects arising in the mesomorphic phase due to colloidal inclusions. Boojums were observed with planar anchoring, whereas Saturn rings were obtained with homeotropic anchoring. The manner in which these topological defects drive colloidal interactions was assessed through a free energy analysis, taking into account the relative orientation between a colloidal dyad and the nematic-field director. The dynamical behavior of defects was qualitatively surveyed from equilibrium trajectories borne from computer simulations.

---

1 Introduction

Nematic liquid crystals with colloidal inclusions provide an avenue for producing multifunctional, tunable materials. The formation of colloidal structures can be controlled through the alignment and colloid-surface anchoring of the host liquid-crystal (LC) solvent.^1–5 Topological defects arise when colloidal inclusions disrupt the prevalent order in the LC solvent. When coupled to the inclusions, these topological defects prompt an elastic binding interaction between guest colloidal particles.^6–10 Defects themselves are key in designing specific colloidal arrangements and in contributing to their stability. Previous studies^11–14 show that these non-covalent binding interactions can be sufficiently strong and capable of withstanding thermal fluctuations, of .

A mounting body of work involving host LC solvents and guest colloidal inclusions has focused on calamitic (i.e., prolate) systems in 2D¹² and 3D.^15,16 The production of sheets,^17,18 wires,^7,15 glasses,¹⁹ crystallites,²⁰ and birefringent soft solids²¹ from bulk dispersions has already been shown. Sensor-based technologies are also possible with colloid–LC composites, for instance, to probe ions,²² point mutations in DNA,²³ and biochemical analytes.^24,25 An alternative approach to the assembly of colloidal particles by LC topological defects involves trapping the inclusions at an interface. Specifically, colloidal structures can be controlled through interfacial disclinations from the LC solvent in contact with an isotropic medium (e.g., air^26–28 or oil²⁹).

Discotic (i.e., oblate) liquid crystals with colloidal inclusions have only been recently reviewed^30,31 though this vein of work began nearly two decades ago. Special consideration has been given to experimental discotic-colloid systems, with an eye to the underlying mesophases and their connection to optical, dielectric, and thermodynamic properties.^32–36 Unlike their calamitic counterparts, discotic systems have not received as much attention in the realm of solvent-assisted colloidal assembly. Research perspectives have focused, instead, on the manner in which nanoparticles enhance overall discotic systems, such as by extending the range of working conditions.^30,31 In fact, discotic systems typically require higher temperatures to achieve the nematic (N_D) phase when compared to the analogous calamitic mesophase.³⁷ However, the need for lower temperatures in practical applications has been previously established.^38–40 A recent experimental breakthrough made it possible to achieve the nematic phase in a discotic system, at room temperature, with gold nanoparticles.⁴¹ A handful of studies have also focused on phase transitions,^42–52 kinetic behavior,⁴⁸ segregation,^53–55 and sedimentation⁵⁶ of colloid-discotic mixtures.

Nematic LC fluids are complex, spatially-oriented systems. Mesogenic domains of LC molecules in an ensemble will align collectively in an average direction known as the director. In the presence of colloidal inclusions, however, an attractive colloid–mesogen interaction can drive an anchoring effect: LC molecules “latch on” to the colloidal surface and undergo a change in orientation with respect to the bulk-field director. Physically, the curvature of colloidal particles coupled with the anisotropic shape of the mesogenic units prevents the LC solvent from occupying volume uniformly near the spatial domain of colloids. The reorientation of LC molecules near the colloid gives rise to fluid regions similar to dislocations, causing defects and elastic distortions.^6–8 As a result, the nematic phase displays a break in symmetry manifested by point or loop defects. Systems can be optimized so that defects “engage” with one another, giving rise to defect-mediated colloidal structures.⁵⁷ The range of possible structures of such LC–colloid composites^58,59 can be expanded with chiral systems.^60–63

Surfactants and thin films^11,64 have been used to condition colloidal surfaces, with the ability to tailor both strength and type of anchoring. In a broad sense, colloidal surfaces can be chemically functionalized to favor either homeotropic (radial) or planar (tangential) anchoring. When planar anchoring is favored, the nematic-field director of the LC medium is expected to be tangential to the colloidal surface, exhibiting a pair of antipodal defects known as boojums that align with the far-field director. On the other hand, when homeotropic anchoring is favored, a loop defect emerges known as a Saturn ring. Another level of control becomes possible with Saturn rings because they can be manipulated with optical tweezers. Additionally, Saturn rings from different colloidal particles can couple to form topological knots,⁶⁵ which obey specific charge conservation laws.¹² In this manner, the observed colloidal structures result from a complex set of interactions prompted by solute–solvent coupling.

Colloidal inclusions cause a (local) distortion with respect to the far-field director, which results in an energy cost taken on by the system. When several inclusions coexist, the system will tend to reduce the volume accessible to these distortions, thus minimizing the cost in solvent elastic energy. This feat is met by restricting colloidal particles to “share” defect zones and makes for the driving force behind colloidal assembly in LC solvents. In calamitic systems, stable aggregates have been obtained with interesting optical features.^66,67 Certain properties have been linked to the resulting symmetries of the assembly, with behaviors analogous to those of photonic crystals^12,68 and with structures consistent with biological systems.^69–71 Our primary goal is to characterize the effective interactions mediated by topological defects in discotic LC solvents with colloidal inclusions.

The present work extends the literature with a Molecular Dynamics (MD) computer simulation study on a discotic system, as a counterpart to the information available for calamitic systems. For example, calamitic systems have been treated along this vein experimentally,⁷² theoretically,^73–76 and with computer simulations.^77–80 An equivalent effort is not currently available, to our knowledge, for discotic systems.² The intent is to draw parallels to and highlight differences on the way topological defects drive colloidal association. In the initial effort reported here, single-colloid and colloid-pair samples are considered as “building blocks” for complex structures.⁸¹

In our model, the “host” solvent is a Gay–Berne discogen, with a parameter set which approximates a triphenylene core. The “guest” colloidal particles are soft spheres. Single-colloid samples are used as a reference system to compare colloid-pair samples. To ensure that attractive interactions in the system are solely due to solute–solvent coupling, and thus to better characterize the emergent interactions, the colloid–colloid potential only takes into account the excluded volume (i.e., it is a soft, purely repulsive interaction). Planar and homeotropic anchoring are compared to a reference arrangement, in which neither type of anchoring is favored. In the spirit of previous studies that considered the relative orientation of colloidal pairs with the far-field director,^82,83 we considered two limiting cases here: when the center-to-center, intercolloidal vector is initialized to be parallel (the PARA case) or perpendicular (the PERP case) to the director. In this manner, four characteristic cases arise: the manner in which topological defects interact is compared between planar (P-PARA and P-PERP) and homeotropic (H-PARA and H-PERP) anchoring. Simulation snapshots and trajectories complement our quantitative results. The treatment of a discotic model, as done in the present work, offers another avenue for structured materials by drawing on the properties of discotic systems.⁸⁴

Our report is structured as follows: in Section 2, the model is reviewed. This includes a brief treatment of the interactions between Gay–Berne discogens, colloidal particles, and the coupling of the two species. Details on the implementation of the model in computer simulations are summarized in Section 3. This section also includes an overview of the computational details for the local order parameter and system free energy. In Section 4, results from computer simulations are presented. Findings are categorized by single-colloid and colloid-pair systems. We qualitatively review the dynamical behavior of topological defects in the four characteristic cases (i.e., P-PARA, P-PERP, H-PARA and H-PERP). The relative stability of the different scenarios is also reviewed in the context of free energy. Closing remarks with a focus on defect-mediated assembly are provided in Section 5.

2 Model

The LC–colloid system contains N_c solid colloidal particles immersed in a nematic LC comprised of N_m mesogens. Position vectors for all colloids are denoted by {R} = R₁,R₂,…,R_Nc, whereas for mesogens making up the solvent, position vectors are denoted by {r} = r₁,r₂,…,r_Nm with their respective (unit) orientation vectors by {ê} = ê₁,ê₂,…,ê_Nm.⁸⁵ The total interaction energy is given by the sum of three contributions,

U_tot({r},{ê},{R}) = U_mm({r},{ê}) + U_mc({r},{ê},{R}) + U_cc({R}), (1)

where U_mm({r},{ê}) is the mesogen–mesogen interaction (dependent on the positions and orientations of all mesogens), U_mc({r},{ê},{R}) is the mesogen–colloid interaction (dependent on the positions and orientations of all mesogens, as well as the positions of all colloidal particles), and U_cc({R}) is the colloid–colloid interaction (dependent only on the positions of all colloidal particles). The last term in eqn (1) is nonzero only when N_c ≥ 2. In our study, we considered LC–colloid systems with N_c = 1 (single-colloid systems) as well as with N_c = 2 (colloid-pair systems). All terms in U_tot({r},{ê},{R}) are detailed in the following subsections.

2.1 Mesogen–mesogen interaction

The mesogen–mesogen contribution U_mm({r},{ê}) assumes pairwise additive interactions. The medium consists of non-spherical, mesogenic particles interacting via the Gay–Berne interaction potential.⁸⁶ An essential feature in tracking thermal effects in the system is the inclusion of both attractive and repulsive interactions, which are captured in the Gay–Berne potential. Each mesogen is represented by an oblate ellipsoid, given that the LC medium is comprised of discotic mesogens (i.e., discogens). The orientation vector for each discogen is defined perpendicular to the plane containing the major axis of the ellipsoid. Because U_mm({r},{ê}) only depends on the relative separation of any two mesogens i and j, the definition is specialized so thatwhere r_ij = r_i − r_j is the center-to-center separation vector between the ith and jth discogen, while ê_i and ê_j are their respective unit orientation vectors. The operational form of the potential is given by

U^mm_ij(r_ij,ê_i,ê_j) = 4ε^mm_ij([Ξ^mm_ij]⁻¹² − [Ξ^mm_ij]⁻⁶). (3)

The relative orientation of mesogens within the medium is taken into account through Ξ^mm_ij and ε^mm_ij, giving rise to anisotropic intermolecular interactions. To define these factors, a fully specified function of a general variable ω is introduced,where c_i = ê_i⋅r̂_ij, c_j = ê_j⋅r̂_ij, c_ij = ê_i⋅ê_j, and r̂_ij = r_ij/|r_ij| is the unit (center-to-center) separation vector. Now,where r_ij = |r_ij| is the magnitude of the separation vector, σ_f is the discogen thickness, and

σ^mm_ij = σ_e[Γ(ω = χ)]^−1/2, (6)

withbeing the length anisotropy of discogens, where κ = σ_f/σ_e is the aspect ratio between the thickness of the discogen and σ_e denotes its diameter. The designations “f” and “e” are mnemonic for face-to-face and edge-to-edge configurations, respectively, as defined in Fig. 1. For discotic mesogens, the aspect ratio is such that σ_f < σ_e and for which we have set σ_e = 1.0.

Fig. 1 Definitions of (a) face-to-face parameters {σ_f, ε_f}, and (b) edge-to-edge parameters {σ_e, ε_e}, associated with the discogen model. Discogens are schematically represented with cylindrical geometry for simplicity, though they are actually oblate ellipsoids in the model.

The energy scale is defined by

ε^mm_ij = ε₀[ε₁(ê_i,ê_j)]^ν[ε₂(r̂_ij,ê_i,ê_j)]^μ, (8)

where ε₀ is the potential energy well depth for two discogens orthogonal to one another (i.e., cross configuration) and to the center-to-center vector (i.e., c_i = c_j = c_ij = 0). Moreover, ν and μ control the contribution of the two dimensionless energy factors, ε₁(ê_i,ê_j) and ε₂(r̂_ij,ê_i,ê_j). In particular,

ε₁(ê_i,ê_j) = [1 − χ²c_ij²]^−1/2, (9)

controls the parallel alignment of discogens and the emergence of mesogenic phases. The second term is given by

ε₂(r̂_ij,ê_i,ê_j) = Γ(ω = χ′), (10)

whereis the energy anisotropy of the discogen, with κ′ = ε_e/ε_f being the ratio of attractive basins of the interaction energy between a face-to-face (ε_f) and an edge-to-edge (ε_e) configuration.

To convey different mesogenic models with the Gay–Berne potential, Bates and Luckhurst⁸⁷ proposed the notation GB(κ, κ′, μ, ν), which is used here: the discogen parameterization in this work is GB(0.345, 0.2, 1.0, 2.0), which corresponds to a (coarse-grained) triphenylene core.^88,89 It has been previously noted^90,91 that κ = 0.345 is higher than triphenylene systems of practical importance due to neglected pendant groups from the discogenic core. Although peripheral groups on the core have been shown in the laboratory to be essential in accessing certain mesophases^84,92–96 (particularly the N_D [discotic nematic] phase), incorporating such detail through a smaller κ in the coarse-grained representation only fine-tunes certain regions of the phase diagram,⁸⁹ without forgoing the essential physics or application features in our studies. The value of κ′ for the current parameterization favors face-to-face over edge-to-edge configurations: this arrangement promotes the nematization of the LC at reasonable thermodynamic state points, representing a key feature in the model for this work. Phase diagrams for GB(0.345, 0.2, 1.0, 2.0) have been previously reported.^97,98

2.2 Mesogen–colloid interaction

The mesogen–colloid interaction is key in the phenomenology of nematic colloid systems. A coupling of length- and timescales effectively occurs, which leads to interesting behaviors behind the break in symmetry of the solvent. The mesogen–colloid interaction depends on the orientation of the interacting mesogens and their relative position with respect to the colloids. In this work, this cross-species interaction is accounted for bywhere X_iα is the separation vector between the ith discogen and the αth colloid (i.e., X_iα = r_i − R_α). Note the separate limits on the summations and the inclusion of the orientation vector ê_i for the ith mesogen.

The mesogen–colloid interaction is captured by a model proposed previously,⁹⁹ which is a generalized potential between spherical and non-spherical particles. This interaction is operationally defined as

within which d_iα represents the separation between the discogen and the colloid surface (i.e., d_iα = X_iα − σ_c/2, where X_iα = |X_iα| is the magnitude of the mesogen–colloid center-to-center separation vector and σ_c denotes the colloid diameter). Anisotropic contributions are captured through σ^mc_iα and ε^mc_iα. Specifically,where

sin(θ) = [1 − cos²(θ)]^1/2 = [1 − c_i²]^1/2, (16)

with c_i defined as in eqn (4), but X̂_iα (the analogous unit vector for r̂_ij) defined using the mesogen–colloid center-to-center separation vector and χ taken from eqn (7). Additionally, σ_e is the discogen diameter defined in the context of eqn (7). Now,where ε_a is the LC–colloid anchoring strength (i.e., an energy scale), μ is matched to that of eqn (8), and

χ′′ = 1 − (κ′′)^1/μ (18)

takes into account the anisotropic contribution in ε_a. In addition, κ′′ = ε_f/ε_e adjusts the ratio of anchoring energy scales on the colloid surface, which effectively controls the anchoring mode: if κ′′ > 1, homeotropic anchoring is favored; if κ′′ < 1, planar anchoring is favored. When κ′′ = 1, planar and homeotropic anchoring are approximately equally favored, as shown in Fig. 2.

Fig. 2 The mathematical form of U^mc_ij(X_iα,ê_i) leads to equal well depths for the LC–colloid limiting configurations shown in the inset only if the colloid diameter . This scenario is relevant when neither planar nor homeotropic anchoring is preferred (i.e., κ′′ = 1). This is confirmed by tracking the minimum of U^mc_ij(X_iα,ê_i) normalized by the energy scale ε_a plotted against σ_c, considering the two limiting anchoring modes [planar (dashed line) and homeotropic (solid line)]. For σ_c = 7σ_e as used in this study, well-depth equivalence for κ′′ = 1 is only approximate (dashed red line). The inset shows U^mc_ij(X_iα,ê_i) as a function of interparticle separation X_iα corresponding to data in the main plot, but specialized to σ_c = 7σ_e and ε_a = 100ε₀. Other anisotropies in anchoring strengths are also shown: κ′′ = 0.5 [red line] and κ′′ = 1.5 [blue line]. Discogens are schematically represented with cylindrical geometry for simplicity, though they are actually oblate ellipsoids in the model. The relative size of colloids and discogens in the schematic are not to scale, but were chosen to emphasize distinct anchoring modes.

To explore the effects of different anchoring strength ratios, we chose κ′′ = {0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 2.0, 5.0, 10.0} for single-colloid systems; the same values were considered for colloid-pair systems with the exception that κ′′ = {4.0, 6.0} was substituted for κ′′ = 5.0 in the set. Additionally, we fixed ε_a = 100ε₀ and σ_c = 7σ_e (i.e., ε₀ = 1.0 and σ_e = 1.0, as indicated in Section 2.1), and μ = 1.0 for all samples.

2.3 Colloid–colloid interaction

The colloid–colloid interaction U_cc({R}) is modeled through a soft-sphere interaction potential, which only accounts for interparticle repulsions on close approach. As such, this term only depends on the separation between colloids and is thus an isotropic contribution defined aswhere R_αβ is the magnitude of the center-to-center separation vector between the αth and βth colloids (i.e., R_αβ = |R_α − R_β|). Furthermore,where ε_c represents the interaction strength between colloid particles, σ_c denotes the colloid diameter, and ϑ controls the softness of the repulsive interaction. For all colloid-pair samples, we used ε_c = ε₀, σ_c = 7σ_e, and ϑ = 36 (i.e., ε₀ = 1.0 and σ_e = 1.0, as indicated in Section 2.1).

As is well known, colloids can undergo clustering or form aggregates through short-range attractive forces (e.g., van der Waals interactions), as well as through dipolar or quadrupolar attractive effects.¹¹ Because our interest is in discerning the ability of the nematic LC in producing effective attractions between colloidal particles, it is important to use the form provided in eqn (20), which neglects any attractions due to colloid–colloid interactions. We thus ensure that any colloid–colloid attractive interactions are solely due to the underlying LC solvent.

3 Methods

Details on implementing the model in computer simulations and the formalisms for analyzing system properties are summarized in this section. As discussed in the opening remarks of Section 2, two types of systems were investigated in this work: single-colloid and colloid-pair samples. In both scenarios, colloidal inclusions were contained in a nematic-phase discotic solvent.

3.1 Molecular dynamics simulations

Trajectories for all samples were acquired with Molecular Dynamics (MD) simulations performed in the canonical (NVT) ensemble. To implement the colloid–LC model, eqn (1) was generalized to include a time variable t: U_tot({r(t)},{ê(t)}, {R(t)}). Simulations were parallelized to expedite the time evolution of all samples. To further optimize the simulations, neighbor lists were used with distance cutoffs as follows: for the mesogen–mesogen contribution, r_cut = 1.6σ_e; for the mesogen–colloid contribution, X_cut = σ_c + σ_f; for the colloid–colloid contribution, R_cut = σ_c + σ_e. Equations of motion were integrated with the velocity-Verlet algorithm, using a simulation time step of δt = 0.001 for translational and orientational dynamics.^100,101 A typical simulation consisted of time steps.

Samples consisted of discogens contained in 3D cubic or rectangular cells. Periodic boundary conditions were applied in all spatial dimensions. Thus, each system represented an unconfined sample of a discotic solvent with colloidal inclusions. For single-colloid samples, N_m = 28 000 discogens were contained in a simulation cell with dimensions L_x = L_y = L_z = 22σ_e. For colloid-pair samples, N_m = 56 000 discogens were contained in a simulation cell with dimensions L_x = L_y = 22σ_e, L_z = 44σ_e. Colloidal dispersions were prepared with a number density ρ = 2.63 and temperature T = 7.0. For this state point, the observed pressure for the bulk solvent is P = 220, located deep in the nematic (N_D) phase region for the Gay–Berne parameterization used in this study.^97,98 The temperature for a given sample was held fixed using a Nosé–Hoover thermostat. A simulation run was performed independently for each value of κ′′ in the study (as described in Section 2.2).

The model presented in this work is in terms of reduced units. To avoid cumbersome notation, all reduced units have been denoted without an asterisk (contrary to convention).¹⁰² Thus, all dimensionalized quantities are reported in “star notation”, e.g., t* refers to a time in ps. As a reference experimental system, we used triphenylene hexa(heptoxybenzoate)¹⁰³ [abbreviated as C7OHBT⁹¹ or H7OBT^94,104] for several reasons: (a) it contains the original mesogenic core (i.e., triphenylene) of the Gay–Berne parameterization,^89,105 (b) it exhibits the discotic nematic (N_D) mesophase relevant in this work,^84,106 and (c) its physically-relevant model parameters have been previously reported. Similar parameter sets have been used to expand the range of studies possible with computer simulations to model specific phenomena in discotic liquid crystals.¹⁰⁷ The coarse-grained representation adopted in this work can be approximately dimensionalized using , , and m* = 1634.064 g mol⁻¹.⁹¹ With these units, key system variables are found as follows: t* = (82.6 ps)t, T* = (194.8 K)T, P* = (1.53 bar)P, and ρ* = (0.0569 nm⁻³)ρ. For example, the dimensionalized simulation time step is approximately equivalent to δt* = (82.6 ps)δt = (82.6 ps)(0.001) = 82.6 fs.

3.2 Orientational order

To characterize interactions that emerge between solvent and colloidal inclusions, the local orientational order of the LC solvent was determined from the Maier–Saupe (nematic) order parameter S_2,loc, obtained by diagonalizing the orientational tensor,where ⊗ denotes the tensor product, I corresponds to the identity matrix, and is the number of mesogens contained in a subvolume of the ensemble [chosen to be sufficiently large such that discogens]. Because our analysis accounts for discogen subpopulations (i.e., through discretized volumes), finite-size corrections are implemented as reported previously.^108–110 For an isotropic sample, S_2,loc = 0. As the population of molecular principal axes (i.e., normal to the plane of the discogen face) aligning with the director in the region n̂_loc increases, S_2,loc will also increase. The extension of eqn (21) to include all N_m discogens in the ensemble yields S₂ = λ_max as the largest eigenvalue (after diagonalizing Q);¹¹¹ the corresponding normalized eigenvector is the director n̂ of the entire sample.

Local orientational order is presented in color maps conveying a transverse plane containing the center-to-center vector of the colloid pair and parallel to one of the faces of the rectangular simulation cell. The local nematic order parameter is averaged over equilibrium simulation steps. Color maps are derived by gridding data from the local orientational order to obtain subregions in the system, with a linear resolution of 0.33σ_e.

3.3 Free energy

For colloid-pair systems, the potential of mean force (PMF) between the colloidal inclusions immersed in a nematic-phase discotic solvent is used to characterize the effective attraction between the dyad. It is stressed that such an attraction is prompted by the underlying solvent only, and not the result of any type of covalent bonding. To perform the calculation, the position vector of a reference colloid (i.e., R₁) is fixed while systematically changing the position vector of the probe colloid (i.e., R₂). The calculation is initiated with a sufficiently large separation R₁₂ = |R₁ − R₂|, such that the force acting on the colloid pair corresponds to that of the bulk solvent (i.e., remains constant). As a next step, R₁₂ is decreased in N_s steps (e.g., this amounts to N_s ∼ 26 steps in our study, from R₁₂ ∼ 9.5σ_e to 7.0σ_e). For each step, the force acting between colloidal inclusions is recorded multiple times from which the mean value is determined at the end of the run. The separation between colloids is constrained for a given value of R₁₂ in order to acquire a statistically meaningful average. The calculation is complete when R₁₂ ≈ σ_c. Integrating the mean force over the trajectory ξ of N_s steps yields the PMF acting between the colloid pair.¹¹²

Because the distance R₁₂ is the order parameter of the PMF, the mean force F(R₁₂) acting between the colloid pair is

where 〈⋯〉_|R1−R2| signifies an average over all configurations in which the separation between R₁ and R₂ matches R₁₂. In addition, the mean value 〈F(R₁₂)〉_|R1−R2| is obtained by averaging over all discogen coordinates and all energy contributions involving the colloid pair. Therefore,where k_B is the Boltzmann constant and the set of differentials dr₁dr₂…dr_Nm in the integral for the canonical average has been shortened to d{r}. Recognizing the denominator of eqn (23) as the configurational partition function Q and that the PMF is related to the Helmholtz free energy A as

A(R₁₂) = −k_BT  ln Q(R₁₂), (24)

it can be established thatAs the distance between colloids R₁₂ is systematically changed (while collecting an average force over the trajectory ξ), each ith step |R₁ − R₂|_i of the N_s steps corresponds to a different system with a different free energy.¹¹² The change in free energy between two steps (e.g., i + 1 and i) is then given by

ΔA(R₁₂)_i+1,i = A(R₁₂)_i+1 − A(R₁₂)_i. (26)

The change in free energy (as a function of R₁₂) is thenwhere dξ represents a trajectory differential [i.e., ξ is the path for the probe colloid (i.e., R₂) to approach the reference colloid (i.e., R₁), so ξ is a function of R₁₂].

4 Results and discussion

Results are presented separately for single-colloid and colloid-pair systems. For colloid-pair systems, the intercolloid (center-to-center) vector R₁₂ can take on an arbitrary orientation with respect to the nematic-field director n̂. Two limiting scenarios are considered here when initializing the system: when R₁₂ is parallel to n̂ (i.e., the PARA case) and when R₁₂ is perpendicular to n̂ (i.e., the PERP case). In both scenarios, R₁₂ is arranged to be parallel to the z-axis (i.e., the longest dimension of the rectangular simulation cell) when showing simulation snapshots, unless otherwise noted. In all cases, the orientation of the far-field director is conveyed in a sphere-referenced coordinate system. The resulting equilibrium arrangement for the two orientation scenarios are treated separately.

4.1 Single-colloid systems

The reference system for single-colloid samples is characterized by discogens that have no preference for either homeotropic or planar anchoring (i.e., κ′′ = 1). Results for the reference system are shown in Fig. 3(a), where no prevalent color is observed for discogens in contact with the colloid surface. Bulk regions of the solvent mostly show the same (green) color, indicating a favored orientation characteristic of the N_D phase. The corresponding orientational order color map in Fig. 3(b) shows that the nematic order parameter fluctuates between S_2,loc = 0.50–0.75, typical of the N_D phase, with no perceptible defect.

Fig. 3 Equilibrium configuration for a single-colloid system with an anchoring strength anisotropy of κ′′ = 1.0. (a) Cross-section in the yz-plane with respect to the colloidal particle in the discotic solvent. Each discogen is colored according to its orientation vector. The average orientation of the far-field director for the system at equilibrium is indicated by the sphere-referenced coordinate system. (b) Color map of the local nematic order parameter S_2,loc, in the same plane as in Panel (a): no topological defect is observed. The color bar is used to characterize S_2,loc within a local region on the color map. Refer to Section 3.2 for computational details of the color map. All colloidal dispersions in this study were prepared with fixed number density (ρ = 2.63) and temperature (T = 7.0), giving a bulk internal pressure of P = 220. Under these conditions, the Gay–Berne model used in this study is located deep in the nematic (N_D) phase.^97,98

Topological defects become pronounced as planar anchoring anisotropy increases (i.e., κ′′ → 0). This can be appreciated in sample cross-sections of Fig. 4 (top panels), comparing data for κ′′ = 0.8 to 0.1. Discogen colors near the colloid surface exhibit planar anchoring (i.e., when discogen edges interact with the colloid surface). When κ′′ = 0.1, the planar anchoring energy overwhelms the competing potential energy of the nematic field: discogens align perpendicular to the director at the equatorial poles of the colloid, giving rise to the formation of boojums. The defect is easily traced in the order parameter color maps of Fig. 4 (bottom panels). The nematic order parameter (S_2,loc ∼ 0.3) is lower for boojums when compared to that of the N_D phase.

Fig. 4 Equilibrium configurations in the yz-plane for an energy anisotropy favoring parallel anchoring: κ′′ = {0.10, 0.20, 0.40, 0.80}. Top panels are as in Fig. 3(a), whereas bottom panels are as in Fig. 3(b). Two opposite regions appear on the colloid surface with orientational disorder (i.e., blue and green shading in the color maps). The effect diminishes as κ′′ → 1. The average orientation of the director in the sphere-referenced coordinate system applies to all samples.

Results for homeotropic anchoring (when κ′′ > 1) are shown in Fig. 5. As expected, sample cross-sections show that discogen faces osculate the colloid surface, an effect that becomes more pronounced as κ′′ increases from unity. The effect of anchoring can be gauged by the color change in discogens near the colloid surface, shown in the top panels of Fig. 5: purple discogen hues are pronounced at the meridian poles of the colloid: the discogens reorient perpendicularly relative to the nematic-field director. Color maps of the nematic order parameter are shown in the bottom panels of Fig. 5. In the bulk regions of the solvent, the orientational order parameter fluctuates around values characteristic of the N_D phase (i.e., S_2,loc ∼ 0.50–0.75). As homeotropic anchoring becomes more favorable, regions with low values of the nematic order parameter (i.e., S_2,loc ∼ 0.25) also increase, indicating the onset of orientational frustration driven by topological defects.

Fig. 5 Same as in Fig. 4, but for an energy anisotropy favoring homeotropic anchoring: κ′′ = {2.0, 5.0, 10.0}. Although the color maps share some similarities to those in Fig. 4, they in fact correspond to Saturn ring defects.

Topological defects from the color maps in Fig. 5 appear to be similar to those in Fig. 4, but are in fact due to a different topological defect: a Saturn ring. Thus, the concentric-like blue shading observed with boojums is lost. An alternate cross-sectional view of the sample is presented in Fig. 6, where the ring structure becomes apparent. Colloidal inclusions immersed in nematic-phase prolate systems (i.e., σ_e < σ_f) for which homeotropic anchoring is favored also display Saturn rings.^79,113 In our systems, topological defects become stable against thermal fluctuations when κ′′ > 2. When 1 ≤ κ′′ ≤ 2, defects are not well defined, as can be inferred from Fig. 6.

Fig. 6 Same as in Fig. 5, but for the xy-plane. Saturn rings are evident in this cross-sectional view.

A visual summary on the topological defects observed in single-colloid systems is provided in Fig. 7. The left panel (i.e., κ′′ < 1) shows that boojums appear when planar anchoring is favored (i.e., κ′′ < 1), which manifest on the colloid surface on opposite ends [cf. Fig. 7(a)]. When there is no preference in anchoring mode (i.e., κ′′ = 1), no topological defects appear in the system [cf. Fig. 7(b)]. Lastly, Saturn rings emerge when homeotropic anchoring is favored (i.e., κ′′ > 1) [cf. Fig. 7(c)].

Fig. 7 A visual summary of topological defects for single-colloid systems in a N_D solvent: (a) boojums form with planar anchoring [i.e., κ′′ < 1], (b) no discernible defects arise when no anchoring type is energetically preferred [i.e., κ′′ = 1], and (c) Saturn rings emerge with homeotropic anchoring [i.e., κ′′ > 1]. The average orientation of the director is shown below each rendition, indicated by the sphere-referenced coordinate system.

4.2 Colloid-pair systems

In this section, colloid-pair samples are treated. All simulation snapshots show an average final (surface-to-surface) distance between colloids of 1.7(±0.2)σ_e. The PARA and PERP cases are treated separately in the following subsections; refer to the opening remarks of Section 4 for a definition of these two cases.

Colloidal dyads do not possess a bond of any sort in this work: each colloid is free to undergo independent translational motion, subject to the excluded volume of other colloidal particles and the discotic solvent. Additionally, a colloidal dyad is not constrained or confined to attain an equilibrium arrangement. Furthermore, there are no external fields used to constrain the relative orientation between the intercolloidal vector and the far-field director attained at equilibrium.

Although the PARA and PERP designations strictly describe initial configurations of different systems, they are still useful when classifying spatial arrangements observed at equilibrium due to the distinctly identifiable behavior in each case. Although the relative orientation between the intercolloidal separation vector and the far-field director (i.e., angle) is conserved, their absolute position in the simulation cell is subject to change. All simulation snapshots shown in this discussion were centered and adjusted so that the intercolloidal vector is parallel to the z-axis of the simulation cell in order to ease comparison between different anchoring modes and arrangements.

4.2.1 Parallel configuration. As was done for single-colloid systems, the reference system in which there is no preference in anchoring (i.e., κ′′ = 1) is considered first. On inspecting the arrangement of discogens throughout the colloid-pair system in Fig. 8(a), there appears to be no preferential ordering on the colloidal surfaces: the solvent occupies the intercolloid region with no preferred orientation. However, the color map of the corresponding configuration, Fig. 8(b), shows that weak but evident defects are present. More specifically, two boojum defects (one set on each colloid) emerge slanted, following the direction of the nematic-field director. Within the intercolloid region, the boojum from one colloid coalesces with that of the other colloid. This behavior is only fleetingly stable when κ′′ = 1, as inferred from a visual inspection of the associated MD trajectory. As for the bulk regions of the system, the color map establishes that discogens, on average, orient along the nematic-field director. The local-field director was assessed for two cases in the limit of strong anchoring: bulk-like homogeneity is recovered beyond the spatial domain of the colloidal-pair inclusions. Sample calculations are deposited as ESI† for this work (see ESI† file Local-Field.pdf).

Fig. 8 As in Fig. 3, but for a colloid-pair immersed in a N_D solvent. The director n̂ is initialized to be parallel to the intercolloid vector R₁₂ (i.e., PARA case). When there is no anchoring preference in the system (i.e., κ′′ = 1), the nematic field stabilizes the colloid dyad through periodic, fleeting coupling of boojum defects.

As planar anchoring becomes more favorable, upon decreasing κ′′ from unity, the colloid surface energy overcomes the field imposed by the nematic phase and the formation of boojums is now stabilized when compared to samples with κ′′ = 1, as was the case for single-colloid systems. We term this arrangement the P-PARA case. An important difference with the colloidal dyad is that the system minimizes the internal energy by having colloids “share” common, orientationally-disordered zones. Configuration snapshots and order parameter color maps for this scenario are shown in Fig. 9. Effectively, the “sharing” of defects leads to attractive interactions between colloidal inclusions in which colloidal assembly is mediated by topological defects.

Fig. 9 As in Fig. 8, but for κ′′ = {0.1, 0.2, 0.4, 0.6} (P-PARA case). The coupling of boojum defects becomes more stable against thermal fluctuations as the anchoring strength anisotropy increasingly favors planar anchoring (i.e., as κ′′ decreases). The average orientation of the director in the sphere-referenced coordinate system applies to all samples.

When increasing κ′′ from unity to induce homeotropic anchoring, an arrangement denoted here as the H-PARA case, Saturn rings appear in a manner consistent with the single-colloid samples. Configuration snapshots and order parameter color maps for this regime are shown in Fig. 10. In similar fashion to the P-PARA case, the Saturn rings in the H-PARA case appear with a slanted geometry, following the orientation of nematic-field director (i.e., a two o'clock slant, versus the ten o'clock slant of the P-PARA case). The difference in how the systems were initialized leads to qualitatively different interactions between defects. In the P-PARA case, defects merge in a sustained manner with time; in the H-PARA case, Saturn rings coalesce in subregions sporadically and transiently.

Fig. 10 As in Fig. 9, but for κ′′ = {2.0, 4.0, 6.0, 10.0} (H-PARA case). Although these snapshots are reminiscent of the coupled boojums in Fig. 8, the defects actually correspond to a pair of slanted Saturn rings. This arrangement is better discerned in Fig. 11.

The slanted arrangement of topological defects, as observed here for the P-PARA or H-PARA cases, has been observed previously in the literature for calamitic systems.^{15,64,83,114–119} Both Saturn rings and boojums are quadrupolar defects¹²⁰ and thus engage with other particles via quadrupolar interactions, obeying characteristic conservation laws.^12,121 This governs the manner in which dimers and higher-order structures form: zigzag or chain-like aggregates are stabilized by exhibiting a slight inclination with respect to the intercolloid vector. In essence, this reduces the volume of the distorted region with respect to the director⁸³ and minimizes the overall energy of the system. A visual summary of the topological defects observed for colloid-pair systems (both P-PARA and H-PARA scenarios) is provided in Fig. 11.

Fig. 11 A visual summary of topological defects for colloid-pair systems initialized with a PARA arrangement, with anchoring anisotropies κ′′ as indicated. Planar anchoring [Panels (a), (b), and (c)] leads to a of coupling boojums. The boojums on each colloid engage even when no specific anchoring is preferred [Panel (d)]. Homeotropic anchoring [Panels (e), (f), and (g)] leads to Saturn ring pairs. The average orientation of the director in the sphere-referenced coordinate system applies to all samples.

To qualitatively assess the temporal behavior of topological defects, we inspected trajectories borne from the MD simulations. Trajectories were visualized and deposited as ESI† for this work. In the P-PARA case, the intercolloid domain is always populated with the topological defect. The intensity with which the defect across the intercolloid domain is shared varies with time. At the periphery of the colloidal dyad, defects are significantly more mobile and their intensity is just as variable (see ESI† file P-PARA.mp4). For the H-PARA case, the shape and slant of the Saturn rings are persistent. With time, the rings fuse transiently along a subregion within the intercolloid domain: this flickering behavior, however, does not yield complete “fusion” of the Saturn ring dyad (see ESI† file H-PARA.mp4).

4.2.2 Perpendicular configuration. In this section, focus is given to the PERP case for colloid-pair samples. Results for both planar and homeotropic anchoring appear in Fig. 12. As expected, the color maps show the appearance of boojums for planar anchoring (i.e., κ′′ = 0.1, 0.4). Specifically, the boojums are mutually parallel and align with the nematic-field director in the P-PERP case. However, the extent of interaction between the set of boojums is sporadic and it is not sustained within the intercolloid gap, as in the P-PARA case. When there is no preference in anchoring (i.e., κ′′ = 1), the system already presents a weak effect that is driven solely by the field of the nematic solvent [cf. Fig. 12(c), as was already seen for the P-PARA case]. For homeotropic anchoring (i.e., κ′′ = 4.0, 10.0), a Saturn ring triad is observed: two coplanar rings perpendicular to n̂ are adjoined by a third ring perpendicular to R₁₂ and to the plane containing the other two rings. The result for the H-PERP case has been previously observed in prolate nematic colloids^15,79 as well as in numerical studies^17,122 and simulations.¹¹³ Shown in Fig. 13 is a visual summary of topological defects for the PERP case for the two anchoring energy anisotropies (i.e., P-PERP and H-PERP).

Fig. 12 As in Fig. 9 and 10, but when the director n̂ is initialized to be perpendicular to the intercolloid vector R₁₂ (PERP case). Planar anchoring [Panels (a) and (b)] leads to a parallel arrangement of boojums (the P-PERP case), which align with the nematic-field director. Even when neither type of anchoring is favored [Panel (c)], the PARA arrangement of boojums is operative. Homeotropic anchoring [Panels (d) and (e)] leads to a triad of Saturn rings: two coplanar rings perpendicular to n̂ adjoined by a third ring perpendicular to both the intercolloid vector R₁₂ and to the plane containing the other two rings. This H-PERP case is better discerned in Fig. 13.

Fig. 13 As in Fig. 11, but for the PERP case. The trends correspond to those already described in Fig. 12. The Saturn ring triad described in Fig. 12 is more evident in Panels (e), (f), and (g). For κ′′ = 1, the formation of boojums is imperceptible in Panel (d), but evident in Panel (c) of Fig. 12: this hints at the highly dynamic nature of the defect.

Similar to the treatment of the dynamical assessment of defects in Section 4.2.1, the P-PERP and H-PERP cases are pursued here. For planar anchoring, the parallel boojum arrangement presents a flickering behavior, but seemingly the intercolloid domain never mediates the topological defect (see ESI† file P-PERP.mp4). Boojums “fuse” transiently from the top or bottom region, in a sideways manner. On the other hand, the Saturn ring triad observed for homeotropic anchoring is stable throughout, showing characteristic solvent fluctuations (see ESI† file H-PERP.mp4). This latter behavior is of interest given that theoretical¹²² and simulation studies with prolate mesogens¹⁷ indicate that the stability of the Saturn ring triad may be metastable or induced by confinement, an effect sensitive to the relative length scales of the system. Although finite-size effects may be operative in these systems, they are still experimentally relevant given that many LC-based applications have spatial design constraints. Additionally, numerical studies have found that the Saturn ring triad is relatively short-lived:¹⁷ we found no indication that this structure is transient on the timescales considered here, nor that it is a time-averaged artifact¹¹³ as shown by the ESI† files.

The dynamical behavior of topological defects, as previously discussed, provide insight into the elastic interactions prompted by the N_D solvent. For instance, fluctuations leading to the transient coalescence of defects across colloidal inclusions is a clear manifestation of the reversible, non-covalent attractive interactions cited as being requisite for molecular self-assembly.^4,5,123 Colloidal structures mediated by LC solvents are thus driven by the system so as to minimize energy penalties. The reversible nature of these interactions is central to self-healing materials,³ in which structural aberrations can dissipate via a fluctuating background imposed by topological defects.

4.3 Free energy of colloid pair inclusions

To characterize the effective interactions present in the colloid-pair samples, the change in free energy ΔA(R₁₂) was determined as a function of intercolloid (center-to-center) separation R₁₂. We remind the reader that the separation between colloids is constrained for a given value of R₁₂ in order to acquire statistically meaningful averages for ΔA(R₁₂). Focus is given to comparing the two limiting geometrical scenarios: the PARA and PERP cases. Results from these calculations are shown in Fig. 14.

Fig. 14 The normalized change in free energy ΔA(R₁₂)ε₀/(k_BT) as a function of the intercolloid (center-to-center) separation R₁₂ for different values of anchoring strength anisotropy: (a) κ′′ = 0.10, (b) κ′′ = 4.0, and (c) κ′′ = 10.0. Shown are PARA (dashed line) and PERP (solid line) cases. As a guide to the eye, the limit ΔA(R₁₂) = 0 is shown (dotted red line). Insets are zoomed regions of ΔA(R₁₂)ε₀/(k_BT) at larger intercolloid distances.

When the colloid dyad engages through planar anchoring (i.e., κ′′ = 0.1), the interaction is slightly more stable [i.e., lower ΔA(R₁₂)] for the PERP case. This can be gauged from the global minimum at R₁₂ ∼ 8σ_e in Fig. 14(a). This length scale corresponds to a surface-to-surface gap between colloids of ∼σ_e, which nominally would accommodate up to a discogen with its edges making contact between the two colloid surfaces. The gap is actually less due to an effective “halo” surrounding the colloid, produced by the mathematical form of the mesogen–colloid interaction (a shell of thickness ∼ 0.3σ_e). Physically, the “halo” is an additional excluded volume allowing only for facial interaction (i.e., the discogen “face” bridges the two colloidal surfaces within the gap), which dominates the global minimum in ΔA(R₁₂). The magnitude of ΔA(R₁₂) for this case is comparable to previously published studies on calamitic samples with colloidal inclusions.⁷⁸ A similar response is observed for the PARA case, albeit with a slightly less stable [i.e., higher ΔA(R₁₂)] global minimum. Sample snapshots for both geometrical scenarios are provided in Fig. 15(a) and (b).

Fig. 15 Configuration snapshots linked to global minima in the free energy profiles of Fig. 14(a) [i.e., when κ′′ = 0.1]. To ease the distinction in the alignment of the director with respect to the intercolloid vector, the color scale for the discogens was adjusted: red samples correspond to the PERP case, whereas green samples are for the PARA case. Panels (a) and (b) are for R₁₂ ∼ σ_e; Panels (c) and (d) are for R₁₂ ∼ 9σ_e. The average orientation of the director in the sphere-referenced coordinate system applies to the two samples in each column.

An additional local minimum is observed (with planar mesogen anchoring) when the intercolloid separation is R₁₂ ∼ 9σ_e: this corresponds to a surface-to-surface gap between colloids of ∼2σ_e. After accounting for the colloidal “halos”, the gap is sufficiently wide to accommodate a discogen positioned so that its edges bridge the gap. Sample snapshots for both scenarios (i.e., PARA and PERP) are shown in Fig. 15(c) and (d). Similar to the response observed for the global minimum, the local minimum is slightly more stable in the PERP case.

The response in ΔA(R₁₂) for 8.2σ_e < R₁₂ < 8.9σ_e arises from the competition between energetically favorable mesogen–colloid interactions and the extent to which orientational order is transmitted from the bulk region of the solvent. In the lower end of the R₁₂ range, there is a slightly higher energy penalty for the PERP case because the discogens in the intercolloid gap are distorted more sharply. This effect reverses in the higher end of the R₁₂ range, where discogens within the gap “blend” their molecular directors with that of the bulk, with a concomitant lowering of ΔA(R₁₂). These scenarios are contrasted in both limiting geometrical scenarios for R₁₂ ∼ 8.4σ_e [Fig. 16(a) and (b)] as well as R₁₂ ∼ 8.8σ_e [Fig. 16(c) and (d)].

Fig. 16 As in Fig. 15, but for the crossover region when 8.2σ_e < R₁₂ < 8.9σ_e. Panels (a) and (b) are for R₁₂ ∼ 8.4σ_e; Panels (c) and (d) are for R₁₂ ∼ 8.8σ_e.

For homeotropic anchoring (i.e., κ′′ = 4.0 and κ′′ = 10), the only prominent feature is a global minimum for both PARA and PERP cases, at a length scale that tends toward intercolloid contact (i.e., R₁₂ ∼ 7σ_e) with increasing κ′′. For the two values of κ′′ considered in this analysis, ΔA(R₁₂) shows enhanced stability in the PARA case. This effect can be traced back to the manner in which Saturn rings are manifested in the two limiting geometries: in the PARA case, two slanted Saturn rings form, wherein transiently a subregion of the rings merges within the intercolloid gap. On the other hand, a Saturn ring triad arises in the PERP case: each colloid is surrounded by its own ring (coplanarly positioned, but perpendicular to the nematic-field director); the two resulting rings are joined by a third one positioned perpendicularly between them. These two scenarios are easily contrasted by revisiting Fig. 11(g) and 13(g). Although the ring-pair structure leads to greater stability [cf. the global minima in Fig. 14(b) and (c)], it has been difficult to compare this structure in the present literature: previous studies have focused on the PERP arrangement. The phenomenology observed in our study for the PERP case is consistent with previous work on calamitic systems.⁷⁸

The study of free energy profiles proves useful in summarizing principles where the anchoring mode becomes a design control. Specifically, if colloidal inclusions are to assemble via the transient coupling of boojums, then an arrangement consistent with PERP geometry enhances the stability of such systems. On the other hand, if the coupling should occur through Saturn rings, the free energy analysis suggests that a PARA arrangement is relatively more stable, even though “communication” between rings is transient. At the expense of some stability, it is possible to attain an adjoined arrangement of Saturn rings with PERP geometry. These comments should be tempered against kinetic studies: it is possible that one arrangement might yield higher-order structures that are more readily attained dynamically. This issue is, however, beyond the scope of the present work.

5 Conclusions

Colloidal inclusions dispersed in a nematic-phase discotic (i.e., oblate mesogen) solvent were studied at the level of Gay–Berne mesogens with Molecular Dynamics simulations. Gay–Berne parameters were tailored to capture a triphenylene core in a coarse-grained manner. Colloids were modeled as soft spheres with a purely repulsive interaction potential. Colloidal particles were coupled to the mesogenic solvent using a generalized interaction potential for mixtures of spherical and nonspherical particles. Systems were investigated containing one or two colloidal inclusions. Equilibrated trajectories from computer simulations were used to characterize the mesomorphic structure in the presence of the colloidal inclusions. Energy strengths in the colloid–mesogen description were adjusted to capture three anchoring modes: planar, homeotropic, and equally favorable.

In all cases, equilibrium configurations were acquired for unconfined systems in the absence of any external fields. The translational motion of colloidal particles was unconstrained when acquiring equilibrium data for solvent-mediated interactions. In the case of colloid-pair samples, each colloidal particle underwent independent motion. Furthermore, the relative arrangement between the intercolloidal separation vector and the far-field director is conserved but subject to typical thermal fluctuations.

Systems possessing a single colloidal inclusion yield either of two topological defects, depending on the anchoring mode. Planar anchoring leads to boojum defects, whereas homeotropic anchoring yields Saturn rings. No discernible topological defect is observed when either of these anchoring modes are equally favored from an energy standpoint, as expected. This is consistent with prior studies focused on calamitic (i.e., prolate mesogen) liquid-crystalline solvents with colloidal inclusions.

A relevant issue when two colloidal inclusions are introduced into the nematic-phase discotic solvent is the effective interaction acting between the colloidal particles, which can lead to self-organization in the host solvent. To characterize the influence of the discotic solvent with respect to the far-field director, two limiting configurations were considered when initializing the system: the intercolloid separation vector of the colloidal pair is either parallel or perpendicular to the director. The pair of colloidal inclusions, when in close proximity, already exhibit boojums even when the colloid–mesogen interaction strengths do not favor one anchoring mode over another. This effect is solely driven by the nematic field of the solvent and is not observed when only one colloidal inclusion is present.

When adjusting the colloid–mesogen interaction strengths to favor either anchoring mode, the effective interaction between colloids will differ depending on the relative alignment of the director. Thus, four scenarios emerge: planar-parallel (P-PARA), planar-perpendicular (P-PERP), homeotropic-parallel (H-PARA), and homeotropic-perpendicular (H-PERP). Boojums are always observed to be axially parallel to the far-field director. Saturn rings appear with homeotropic anchoring and couple in a manner dependent on the relative arrangement of the colloid pair. Saturn rings arrange in parallel planes perpendicular to the far-field director. A triad of Saturn rings is observed in the H-PERP case: two coplanar rings perpendicular to the director are adjoined by a third ring that forms in the intercolloid domain, but parallel to the director.

The manner in which the nematic-phase discotic solvent drives an effective interaction between colloidal particles was assessed with free energy profiles. The order parameter for this calculation is the intercolloidal distance: the relative distance between a dyad is fixed for each value in the sweep of the order parameter to obtain statistically sound averages. When boojum defects arise in the system, colloids interact more favorably in the P-PERP case. When Saturn rings are involved, the interaction between colloids is more favorable in the H-PARA case. This orientation has been difficult to trace in the literature since most studies have previously focused on perpendicular arrangements.

Configuration snapshots and free energy profiles used as aids to characterize mesomorphic behavior are static in nature. Equally important is having access to dynamical information of such systems. To glean qualitative dynamical information of the topological defects due to colloidal-pair inclusions, trajectories from Molecular Dynamics simulations were inspected. Temporal windows were visualized from equilibrated samples for the four scenarios previously discussed. With the exception of the Saturn ring triad observed in the H-PERP case, defects display high spatial and temporal variation, and thus appear to “flicker”. In the P-PARA case, boojums merge in the intercolloid region. In the P-PERP case, the intercolloid region is not populated; instead, boojums merge periodically, for brief episodes, in a side-to-side fashion. The Saturn ring dyad in the H-PARA case, in which each colloid possesses its own ring, merges over a small region within the intercolloidal gap, briefly and periodically. The Saturn ring triad appearing in the H-PERP case yields a structure that persists spatially, with only thermal fluctuations being evident.

The ability to control topological defects is a relevant design parameter in many liquid-crystalline-based technologies. As our results show, the extent to which topological defects engage is a sensitive function of both anchoring mode and the relative arrangement of the nematic-field director with respect to colloidal inclusions. This work has combined static and dynamic information to better establish how effective interactions prompted by the nematic solvent aid in structuring colloidal inclusions. Our focus has been on a discotic solvent, given that it has not been previously used in the context of the issues we pursued here. The versatility of discotic systems, however, can extend the range of application given their unique mesophases and wider working range. It is anticipated that new modes to colloidal association may be possible in discotic solvents, in part because the molecular aspect ratio imparts a phase morphology different from that of calamitic systems. Future work in our group will focus on extending the complexity of colloidal assembly and quantitative dynamical analysis attainable from computer simulations.

6 Description of the ESI

The orientation of the local-field director throughout a sample was explored in the limit of strong anchoring for the PARA geometry. Results are reported in ESI† file Local-Field.pdf. Other anchoring strengths as well as the PERP cases explored in this work yield comparable results. The nematic field recovers bulk-like homogeneity beyond the spatial domain of colloidal inclusions.

Trajectories from equilibrated samples were visualized over a sufficiently long time interval to capture characteristic fluctuations of topological defects in nematic LC–colloid systems, for four characteristic scenarios. In each scenario, only the topological defect and the colloidal inclusions are shown for clarity; views were chosen to highlight the manner in which defects engage with a colloidal pair. Each video is differentiated by its anchoring mode and relative orientation (i.e., between the nematic field director n̂ and the intercolloidal separation vector R₁₂):

• H-PARA.mp4.

Homeotropic anchoring (κ′′ = 10.0), when n̂ is initialized to be parallel to R₁₂.

• H-PERP.mp4.

Homeotropic anchoring (κ′′ = 10.0), when n̂ is initialized to be perpendicular to R₁₂.

• P-PARA.mp4.

Planar anchoring (κ′′ = 0.10), when n̂ is initialized to be parallel to R₁₂.

• P-PERP.mp4.

Planar anchoring (κ′′ = 0.10), when n̂ is initialized to be perpendicular to R₁₂.

Although colloidal inclusions do not appear to undergo significant displacements in the videos, each colloid experiences typical thermal fluctuations imposed by the underlying discotic solvent. Translational motion of each colloid is entirely unconstrained (i.e., there are no covalent bonds between colloids that would keep their relative distance and/or orientation with respect to the far-field director constant).

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge the support granted by Red Temática de Materia Condensada Blanda, under the auspices of CONACYT-México. A. D. G.-M. is grateful for a doctoral research fellowship funded through Universidad Autónoma Metropolitana (UAM). J. A. M.-R. acknowledges financial support received from “Fronteras de la Ciencia” CONACYT-México (Grant No. 2015-02-1450). Computational resources were provided by Universidad Nacional Autónoma de México (Grant No. LANCAD-UNAM-DGTIC-276).

Notes and references

1. P. van der Asdonk and P. H. J. Kouwer, Chem. Soc. Rev., 2017, 46, 5935–5949.
2. S. H. Eichhorn and J. K. Yu, in Directed Assembly and Self-organization of Metal Nanoparticles in Two and Three Dimensions, ed. Q. Li, Springer International Publishing, Cham, 2015, pp. 289–336.
3. O. Stamatoiu, J. Mirzaei, X. Feng and H. Torsten, in Nanoparticles in Liquid Crystals and Liquid Crystalline Nanoparticles, ed. C. Tschierske, Springer, Berlin, Heidelberg, 2012, pp. 331–393.
4. K. Thorkelsson, P. Bai and T. Xu, Nano Today, 2015, 10, 48–66.
5. X. Wang, D. S. Miller, E. Bukusoglu, J. J. de Pablo and N. L. Abbott, Nat. Mater., 2016, 15, 106–112.
6. F. Brochard and P. G. de Gennes, J. Phys. (France), 1970, 31, 691–708.
7. P. Poulin, H. Stark, T. C. Lubensky and D. A. Weitz, Science, 1997, 275, 1770–1773.
8. T. C. Lubensky, D. Pettey, N. Currier and H. Stark, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 610–625.
9. M. Kleman and O. D. Lavrentovich, Philos. Mag., 2006, 89, 4117–4137.
10. I. Muševič, Materials (Basel), 2018, 11, 24.
11. M. Tasinkevych, N. M. Silvestre and M. M. Telo da Gama, New J. Phys., 2012, 14, 073030.
12. I. Muševič, M. Škarabot, U. Tkalec, M. Ravnik and S. Žumer, Science, 2006, 313, 954–958.
13. J.-C. Loudet, P. Barois and P. Poulin, Nature, 2000, 407, 611–613.
14. A. Nych, U. Ognysta, M. Škarabot, M. Ravnik, S. Žumer and I. Muševič, Nat. Commun., 2013, 4, 1489.
15. M. Ravnik, M. Škarabot, S. Žumer, U. Tkalec, I. Poberaj, D. Babič, N. Osterman and I. Muševič, Phys. Rev. Lett., 2007, 99, 247801.
16. M. Ravnik, G. P. Alexander, J. M. Yeomans and S. Žumer, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 5188–5192.
17. T. Araki and H. Tanaka, Phys. Rev. Lett., 2006, 97, 127801.
18. N. Hijnen, T. A. Wood, D. Wilson and P. S. Clegg, Langmuir, 2010, 26, 13502–13510.
19. T. A. Wood, J. S. Lintuvuori, A. B. Schofield, D. Marenduzzo and W. C. K. Poon, Science, 2011, 334, 79–83.
20. N. M. Silvestre, Q. Liu, B. Senyuk, I. I. Smalyukh and M. Tasinkevych, Phys. Rev. Lett., 2014, 112, 225501.
21. S. P. Meeker, W. C. K. Poon, J. Crain and E. M. Terentjev, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 61, R6083.
22. Y. Garbovskiy and I. Glushchenko, Appl. Phys. Lett., 2015, 107, 041106.
23. A. D. Price and D. K. Schwartz, J. Am. Chem. Soc., 2008, 130, 8188–8194.
24. J. G. M. Koenig, I.-H. Lin and N. L. Abbott, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 3998–4003.
25. I.-H. Lin, D. S. Miller, P. J. Bertics, C. J. M. J. J. de Pablo and N. L. Abbott, Science, 2011, 332, 1297.
26. A. B. Nych, U. M. Ognysta, V. M. Pergamenshchik, V. G. N. I. Muševič, M. Škarabot and O. D. Lavrentovich, Phys. Rev. Lett., 2007, 98, 057801.
27. M. A. Gharbi, M. Nobili, M. In, G. Prévot, P. Galatola, J.-B. Fournier and C. Blanc, Soft Matter, 2011, 7, 1467–1471.
28. J. M. Cavallaro, M. A. Gharbi, D. A. Beller, S. Čopar, Z. Shi, T. Baumgart, S. Yang, R. D. Kamien and K. J. Stebe, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 18804–18808.
29. J. S. Lintuvuori, A. C. Pawsey, K. Stratford, M. E. Cates, P. S. Clegg and D. Marenduzzo, Phys. Rev. Lett., 2013, 110, 187801.
30. S. Kumar, NPG Asia Mater., 2014, 6, e82.
31. A. Gowda and S. Kumar, Materials (Basel), 2018, 11, 382.
32. P. Yaduvanshi, A. Mishra, S. Kumar and R. Dhar, J. Mol. Liq., 2015, 208, 160–164.
33. M. Mishra, S. Kumar and R. Dhar, RSC Adv., 2014, 4, 62404–62412.
34. S. Varshney, M. Kumar, A. Gowda and S. Kumar, J. Mol. Liq., 2017, 238, 290–295.
35. M. Kumar, S. Varshney, A. Gowda and S. Kumar, J. Mol. Liq., 2017, 241, 666–674.
36. M. Mishra, Phase Transitions, 2018, 91, 752–758.
37. H. K. Bisoyi and S. Kumar, Chem. Soc. Rev., 2010, 39, 264–285.
38. M. Gupta, S. P. Gupta, M. V. Rasna, D. Adhikari, S. Dharaand and S. K. Pal, Chem. Commun., 2017, 53, 3014–3017.
39. H.-H. Chen, H.-A. Lin, S.-C. Chien, T.-H. Wang, H.-F. Hsu, T.-L. Shih and C. Wu, J. Mater. Chem., 2012, 22, 12718–12722.
40. S. Kohmoto, E. Mori and K. Kishikawa, J. Am. Chem. Soc., 2007, 129, 13364–13365.
41. M. Gupta, S. S. Mohapatra, S. Dhara and S. K. Pal, J. Mater. Chem. C, 2018, 6, 2303–2310.
42. S. K. Lai and T. T. T. Van, Colloids Surf., A, 2017, 532, 270–281.
43. M. Chen, H. Li, Y. C. A. F. Mejia, X. Wang and Z. Cheng, Soft Matter, 2015, 11, 5775–5779.
44. D. Kleshchanok, A. V. Petukhov, P. Holmqvist, D. V. Byelov and H. N. W. Lekkerkerker, Langmuir, 2010, 26, 13614–13621.
45. Á. González García, R. Tuinier, J. V. Maring, J. Opdam, H. H. Wensink and H. N. W. Lekkerkerker, Mol. Phys., 2018, 116, 2757–2772.
46. F. Gámez, R. D. Acemel and A. Cuetos, Mol. Phys., 2013, 111, 3136–3146.
47. M. Chen, M. He, P. Lin, Y. Chen and Z. Cheng, Soft Matter, 2017, 13, 4457–4463.
48. T. Schilling, S. Dorosz, M. Radu, M. Mathew, S. Jungblut and K. Binder, Eur. Phys. J.: Spec. Top., 2013, 222, 3039–3052.
49. L. Harnau and S. Dietrich, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 051501.
50. N. Doshi, G. Cinacchi, J. S. van Duijneveldt, T. Cosgrove, S. W. Prescott, I. Grillo, J. Phipps and D. I. Gittins, J. Phys.: Condens. Matter, 2011, 23, 194109.
51. S. M. Oversteegen and H. N. W. Lekkerkerker, J. Chem. Phys., 2004, 120, 2470–2474.
52. R. Aliabadi, M. Moradi and S. Varga, J. Chem. Phys., 2016, 144, 074902.
53. D. de las Heras and M. Schmidt, Philos. Trans. R. Soc., A, 2013, 371, 20120259.
54. J. Sui, Soft Matter, 2019, 15, 4714–4722.
55. D. Kleshchanok, J.-M. Meijer, A. V. Petukhov, G. Portale and H. N. W. Lekkerkerker, Soft Matter, 2012, 8, 191–197.
56. D. de las Heras, N. Doshi, T. Cosgrove, J. Phipps, D. I. Gittins, J. S. van Duijneveldt and M. Schmidt, Sci. Rep., 2012, 2, 789.
57. R. P. Trivedi, I. I. Klevets, B. Senyuk, T. Lee and I. I. Smalyukh, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 4744–4749.
58. R. Bitar, G. Agez and M. Mitov, Soft Matter, 2011, 7, 8198–8206.
59. M. Mitov, C. Portet, C. Bourgerette, E. Snoeck and M. Verelst, Nat. Mater., 2002, 1, 229–231.
60. Y. Yuan, A. Martinez, B. Senyuk, M. Tasinkevych and I. I. Smalyukh, Nat. Mater., 2018, 17, 71–79.
61. I. Muševič, Eur. Phys. J.: Spec. Top., 2019, 227, 2455–2485.
62. K. Nayani, Y.-K. Kim and N. L. Abbott, Nat. Mater., 2018, 17, 14–15.
63. I. I. Smalyukh, Annu. Rev. Condens. Matter Phys., 2018, 9, 207–226.
64. M. Škarabot, M. Ravnik, S. Žumer, U. Tkalec, I. Poberaj, D. Babič, N. Osterman and I. Muševič, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 031705.
65. U. Tkalec and I. Muševič, Soft Matter, 2013, 9, 8140–8150.
66. J. Müller, C. Sönnichsen, H. von Poschinger, G. von Plessen, T. A. Klar and J. Feldmann, Appl. Phys. Lett., 2002, 81, 171–173.
67. H. Qi and T. Hegmann, J. Mater. Chem., 2008, 18, 3288–3294.
68. J. D. Joannopoulos, P. R. Villeneuve and S. Fan, Nature, 1997, 386, 143–149.
69. J. Teyssier, S. V. Saenko, D. van der Marel and M. C. Milinkovitch, Nat. Commun., 2015, 6, 6368.
70. P. Vukusic and J. R. Sambles, Nature, 2003, 424, 852–855.
71. V. Saranathan, C. O. Osuji, S. G. J. Mochrie, H. Noh, S. Narayanan, A. Sandy, E. R. Dufresne and R. O. Prum, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 11676–11681.
72. K. Lin, J. C. Crocker, A. C. Zeri and A. G. Yodh, Phys. Rev. Lett., 2001, 87, 088301.
73. P. Galatola and J.-B. Fournier, Phys. Rev. Lett., 2001, 86, 3915–3918.
74. P. Galatola, J.-B. Fournier and H. Stark, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 031404.
75. E. M. Terentjev, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 51, 1330–1337.
76. B. I. Lev, S. B. Chernyshuk, P. M. Tomchuk and H. Yokoyama, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 021709.
77. D. Andrienko, M. Tasinkevych, P. Patrício, M. P. Allen and M. M. Telo da Gama, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 051702.
78. E. Kim, O. Guzmán, S. Grollau, N. L. Abbott and J. J. de Pablo, J. Chem. Phys., 2004, 121, 1949–1961.
79. O. Guzmán, E. B. Kim, S. Grollau, N. L. Abbott and J. J. de Pablo, Phys. Rev. Lett., 2003, 91, 235507.
80. M. S. Al-Barwani, G. S. Sutcliffe and M. P. Allen, J. Phys. Chem. B, 2004, 108, 6663–6666.
81. Z. Dogic, P. Sharma and M. J. Zakhary, Annu. Rev. Condens. Matter Phys., 2014, 5, 137–157.
82. R. W. Ruhwandl and E. M. Terentjev, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 55, 2958–2961.
83. M. R. Mozaffari, M. Babadi, J. Fukuda and M. R. Ejtehadi, Soft Matter, 2011, 7, 1107–1113.
84. A. N. Cammidge and R. J. Bushby, in Handbook of Liquid Crystals, Part II: Discotic Liquid Crystals, ed. D. Demus, J. Goodby, G. W. Gray, H.-W. Spiess and V. Vill, Oxford University Press, Weinheim, 1998, vol. 2B, ch. VII: Synthesis and Structural Features, pp. 703–705.
85. All unit vectors in this work are denoted with “hat notation”, e.g., r̂_ij = r_ij/|r_ij| is the unit (center-to-center) separation vector between the ith and jth particles.
86. J. G. Gay and B. J. Berne, J. Chem. Phys., 1981, 74, 3316–3319.
87. M. A. Bates and G. R. Luckhurst, J. Chem. Phys., 1999, 110, 7087–7108.
88. A. P. J. Emerson, G. R. Luckhurst and S. G. Whatling, Mol. Phys., 1994, 82, 113–124.
89. D. Caprion, L. Bellier-Castella and J. P. Ryckaert, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 041703.
90. S. Orlandi, L. Muccioli, M. Ricci, R. Berardi and C. Zannoni, Chem. Cent. J., 2007, 1, 15.
91. R. H. Hurt and Y. Hing, Carbon, 1999, 37, 281–292.
92. N. H. Tinh, H. Gasparoux and C. Destrade, Mol. Cryst. Liq. Cryst., 1981, 68, 101–111.
93. P. Hindmarsh, M. Hird, P. Styring and J. W. Goodby, J. Mater. Chem., 1993, 3, 1117–1128.
94. T. J. Phillips, J. C. Jones and D. G. McDonnell, Liq. Cryst., 1993, 15, 203–215.
95. P. Hindmarsh, M. J. Watson, M. Hird and J. W. Goodby, J. Mater. Chem., 1995, 5, 2111–2123.
96. R. J. Bushby, K. J. Donovan, T. Kreouzis and O. R. Lozman, Opto-Electron. Rev., 2005, 13, 269–279.
97. O. Cienega-Cacerez, J. A. Moreno-Razo, E. Díaz-Herrera and E. J. Sambriski, Soft Matter, 2014, 10, 3171–3182.
98. 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.
99. D. Antypov and D. J. Cleaver, J. Phys.: Condens. Matter, 2004, 16, S1887–S1900.
100. J. M. Ilnytskyi and M. R. Wilson, Comput. Phys. Commun., 2001, 134, 23–32.
101. J. M. Ilnytskyi and M. R. Wilson, Comput. Phys. Commun., 2002, 148, 43–58.
102. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 2nd edn, 2017, pp. 487–488.
103. The system used to dimensionalize the Gay–Berne model in this work can be found in chemical repositories as hexakis(4-heptyloxybenzoic acid)triphenylene-2,3,6,7,10,11-hexyl ester [InChI Key: MQEIKEUOALPXJH-UHFFFAOYSA-N. Formula: C₁₀₂H₁₂₀O₁₈].
104. N. H. Tinh, C. Destrade and H. Gasparoux, Phys. Lett. A, 1979, 72, 251–254.
105. M. A. Bates and G. R. Luckhurst, J. Chem. Phys., 1996, 104, 6696–6709.
106. Several experimental systems based on hexaalkoxy benzoates of triphenylene⁸⁴ show that the discotic nematic (N_D) phase is attainable.
107. I. R. Thompson, M. K. Coe, A. B. Walker, M. Ricci, O. M. Roscioni and C. Zannoni, Phys. Rev. Mater., 2018, 2, 064601.
108. R. Eppenga and D. Frenkel, Mol. Phys., 1984, 52, 1303–1334.
109. A. Richter and T. Gruhn, J. Chem. Phys., 2006, 125, 064908.
110. G. D. Wall and D. J. Cleaver, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 56, 4306.
111. Local orientational order (as defined by the Maier–Saupe order parameter) is accounted for by slabs parallel to the xy-plane, traced along the z-axis. An analogous definition can be extended to other spatial dimensions, as needed.
112. W. Billes, F. Bazant-Hegemark, M. Mecke, M. Wendland and J. Fischer, Langmuir, 2003, 19, 10862–10868.
113. A. Humpert, S. F. Brown and M. P. Allen, Liq. Cryst., 2018, 45, 59–69.
114. P. Poulin and D. A. Weitz, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 626.
115. I. I. Smalyukh, O. D. Lavrentovich, A. N. Kuzmin, A. V. Kachynski and P. N. Prasad, Phys. Rev. Lett., 2005, 95, 157801.
116. J. Kotar, M. Vilfan, N. Osterman, D. Babic, M. Copic and I. Poberaj, Phys. Rev. Lett., 2006, 96, 207801.
117. M. Ravnik and S. Žumer, Liq. Cryst., 2009, 36, 1201–1214.
118. M. Škarabot and I. Muševič, Soft Matter, 2010, 6, 5476–5481.
119. M. Škarabot, A. V. Ryzhkova and I. Muševič, J. Mol. Liq., 2018, 267, 384–389.
120. C. Blanca, D. Coursault and E. Lacaze, Liq. Cryst. Rev., 2013, 1, 83–109.
121. H. Stark, Phys. Rep., 2001, 351, 387–474.
122. Y. Wang, P. Zhang and J. Z. Y. Chen, Phys. Rev. E, 2017, 96, 042702.
123. G. M. Whitesides and B. Grzybowski, Science, 2002, 295, 2418–2421.

---

Footnote

† Electronic supplementary information (ESI) available: Additional simulation results and molecular trajectories for characteristic cases are available; refer to the “Description of the ESI” section for a full description. See DOI: 10.1039/c9ra05377h

---