Self-Assembly of an Equimolar Mixture of Liquid Crystals and Magnetic Nanoparticles

: We studied the equilibrium self-assembly of an equimolar mixture of uniaxial liquid crystals (LCs) and magnetic nanoparticles (MNPs) using molecular dynamics simulations. The LCs are modeled by ellipsoids interacting via Gay–Berne potential, and MNPs are represented by dipolar soft spheres (DSS). We found that the LCs show isotropic, nematic, and smectic phases when the mixture is compressed at a ﬁxed temperature. The DSS form chain-like structures, which remain randomly oriented at low densities where the LCs are in the isotropic phase. At intermediate and high densities, the DSS chains align along the nematic and smectic directors of LCs. We found that the DSS inside a chain follow a ferromagnetic ordering. However, the mixture does not show a signiﬁcant macroscopic magnetization. The extent of nematic order in the DSS remains very similar to the LCs in intermediate densities. At high densities, the DSS have a lower extent of nematic order than the LCs. The structure of the LC–DSS mixture was further analyzed via projected pair correlation functions for distances parallel and perpendicular to directors in the nematic and smectic phases.


Introduction
Inclusions of nanoparticles in a liquid crystalline matrix have attracted much research attention in the past few years [1,2]. These novel hybrid soft materials find interesting technical applications in electro-optical devices, material design, medicine, and biotechnology [3][4][5]. In 1970, Brochard and de Gennes proposed that adding a small number of magnetic nanoparticles (MNPs) can facilitate nematic ordering in liquid crystals (LCs) [6]. Since then, several experimental studies on the mixtures of lyotropic and thermotropic liquid crystals of various shapes and magnetic nanoparticles have been reported elucidating the magnetic-field-induced nematic order in LCs [7][8][9][10]. The primary focus of these early studies was on the enhanced nematic ordering in the LC matrix due to the inclusion of magnetic nanoparticles. In contrast, much of the recent studies have focused on using MNPs to introduce functionalities to liquid crystals [5], as the coupling of magnetic nanoparticles to the external magnetic field [11] can be utilized to control the orientation of the LC director.
The self-assembly of the LC-MNP mixtures and their response to external fields have been well explored in experiments [12][13][14][15][16][17][18][19][20][21][22][23]. Recent experiments have demonstrated that the mixture of ferromagnetic platelets suspended in nematic LCs can have stable magnetic domains, which makes these mixtures useful for potential magneto-optical devices due to the switching ability of magnetic domains under the external magnetic field [22,24]. Furthermore, the mechanism of domain growth in these ferronematic mixtures is mediated via the magneto-nematic coupling between MNPs and LCs [25]. Using Landau-de Gennes's free energy formalism, recent theoretical studies have shown that mixtures with symmetric magneto-nematic coupling display a slow growth of domains. In contrast, subdomain morphologies dominated by interfacial defects can occur for asymmetric coupling strengths [26,27].
Recent computer simulations and experimental studies have demonstrated that the small inclusion of magnetic nanoparticles strongly affects the isotropic to nematic (I-N) phase transition in LCs in the presence of an external magnetic field [34,38,39]. The mixtures with large volume fractions of MNPs are expected to display more intriguing ordering as soft-sphere repulsive magnetic nanoparticles display a spontaneous magnetization at high densities. The self-assembly of such mixtures has rarely been explored in computer simulations and is absolutely necessary in view of the fact that the concentration of the nanoparticles strongly affects the I-N transition in these hybrid mixtures [40,41].
In the present work, we considered an equimolar mixture of uniaxial LCs and MNPs of similar size and analyze the self-assembly of the LCs and MNPs as the mixture is compressed at a fixed temperature. We modeled uniaxial LCs by ellipsoids interacting via the generalized Gay-Berne potential [42][43][44][45][46] and MNPs by repulsive spheres with embedded point dipole, the so-called dipolar soft spheres (DSS) [47]. In equilibrium, the DSS form chain-like structures, which remain randomly oriented at low densities while showing spontaneous magnetization at high densities [47][48][49]. We found that the LC-DSS mixture shows the isotropic phase at low densities where the LCs and DSS chains remain randomly orientated. At intermediate and high densities, the LCs display the nematic and smectic B phases, respectively. The DSS follow the orientation of the LCs in the nematic and smectic phases. The overall magnetization of the LC-DSS mixture remains very low, although the volume fraction of the DSS is large at high densities. The structure of the LC-DSS mixture was further analyzed in different phases using pair correlation functions.
The rest of the paper is organized as follows. The details of the model and simulation are given in Section 2. In Section 3, we present various order parameters characterizing the orientational and ferromagnetic ordering in the LC-MNP mixture. We further discuss the pair correlation functions for each species of the mixtures in different phases. Finally, in Section 4, we summarize our findings and discuss the future outlook of the work.

Simulation Details
We considered a binary mixture of LC and DSS with a composition ratio of 50:50 and performed MD simulations in the NVT ensemble using the LAMMPS package [44,50]. The details of the model are the same as used in our previous works [36,51]. The LCs were modeled by ellipsoids interacting via a generalized Gay-Berne (GB) potential, which is defined as follows [44][45][46]: where A i represents the rotation matrix for a particle i. Further, r ij is the center-to-center distance vector between particles i and j, andr ij is the unit vector along r ij . The function U r A i , A j , r ij controls the distance dependence of the GB potential, and the second and third terms in Equation (1) define the dependence of the interaction potential on the relative orientation and distance of the ellipsoids. We considered uniaxial LCs with the lengths of the semi-axes, σ e a = σ e b = σ 0 and σ e c = 3σ 0 ; see Figure 1a. These length parameters set the aspect ratio of LCs equal to 3. The energy parameters for ellipsoids were chosen such that the energy of side-to-side configuration would be five times stronger than end-to-end configuration. The units of the length and energy were set by the width of the ellipsoid, σ 0 , and the well-depth of the GB potential, 0 . The cutoff radius was set to r GB c = 4.0σ 0 , and the empirical parameters of GB potential were set to µ = 1.0, ν = 2.0 and γ = 1.0. The parameter µ is the same as µ used for GB in Brown et al. [44]. Here, reserve µ to denote the dipole moment of the DSS.
The interaction among the LC and DSS was also modeled by a GB potential, with shape and energy parameters for DSS taken as σ s a = σ s b = σ s c = σ 0 and s a = s b = s c = 0 . The cutoff radius was set the same as r GB c . σ 0 , defines the unit of the length. In order to set the aspect ratio equal to 3, the end-to-end length of an ellipsoid is taken to be equal to 3σ 0 . Panel (b): the schematic diagram showing two dipolar soft spheres with embedded permanent dipole moment represented by vectors µ i and µ j . The vector r ij denotes the center-to-center distance between the two dipolar soft spheres. The diameter of the DSS is considered to be equal to σ 0 , the width of the ellipsoids.
The DSS of diameter σ 0 (see Figure 1b) interact via a combination of a soft sphere potential and dipolar interactions [47,48]. The full potential for two DSS particles i and j with dipole moments µ i and µ j is defined as [48] where U SR r ij is the shifted-force soft sphere interaction given as in Equation (3), the cutoff radius for soft sphere potential is set to r SS c = 2.5σ 0 , and U SS r ij = 4 0 σ 0 /r ij 12 .
The long range dipolar interactions are treated with the three dimensional Ewald sum [49,52]. The parameters that characterize the structure and phase behavior of the mixture are the reduced temperature T * = k B T/ 0 , the reduced number density ρ * = Nσ 3 0 /V (where N and V are the total number of particles and total volume, respectively), and the reduced dipole moment µ * = µ 2 / 0 σ 3 . Newton's equations of motion for force and torque were integrated into the NVT ensemble using the velocity Verlet algorithm. A reduced MD time step ∆t * = ∆t/ mσ 2 0 / 0 = 0.002 was chosen, and simulations were performed at fixed µ * = 3.0 and T * , and various ρ * . The value of the dipolar coupling parameter λ = µ 2 /k B Tσ 3 0 is 9.0. At this value of λ, the DSS form long chain-like structures [48,53]. We considered a system consisting of 1000 LC (N e ) and 1000 DSS particles (N s ). We started with a mixture equilibrated at low density and high temperature and then quenched it to the desired T * . Subsequently, we slowly compressed the mixture to obtain desired densities, keeping T * constant by applying a Langevin thermostat [54]. At each compression step, the mixture was equilibrated for 4 × 10 6 time steps. All results presented were averaged over 20 independently prepared samples.

Results
We considered the LC-DSS mixture at a low density, ρ * = 0.1, and temperature, T * = 1.0, and slowly compressed it to achieve different densities in the range. To identify different orientationally ordered phases, we investigated the variation of the nematic ordering of each component of the LC-DSS mixture as a function of density. The nematic ordering was quantified by the nematic parameter, which is defined as the largest eigenvalue of the ordering tensor, namely, the Q-tensor. For the LCs, the components of the Q-tensor are defined as [54,55]: where α, β = x, y, z, and u i α are the components of the unit vector,û i , representing the orientation of the ith ellipsoid. For the DSS, the components of the Q-tensor are defined as [47,48] where µ i α are the components of the unit dipole vector,μ i , for the ith DSS. Note that due to the heat-to-tail invariance of the ellipsoids,û is equivalent to −û, in contrast to the DSS where µ = −µ.
The nematic order parameters S e and S s for the LCs and DSS are defined as the largest eigenvalue of the corresponding Q matrix defined in Equations (5) and (6) [55,56]. The normalized eigenvectors corresponding to S e and S s define the directorsn e andn s for LC and DSS, respectively. We note that the directors are unit vectors representing the average alignment of the LCs and DSS andn e,s is physically equivalent to −n e,s . The isotropic phase of a liquid crystalline system is characterized by a low value (<0.4) of the nematic order parameter, while the I-N transition occurs when the order parameter attains a value of ∼0.4 [56]. Figure 2 shows the variation of the nematic order parameters S e and S s , for the LC and DSS, respectively, as a function of density at temperature T * = 1.0. At low densities, both S e and S s remain close to zero, start to increase around ρ * = 0.4, and attain high values close to 1 at higher densities. This behavior is similar to LC-DSS mixtures with an 80:20 ratio [33], except for the fact that the I-N transition in the present case is observed at higher densities for T * = 1. The nematic phase lies in the range ρ * = 0.4 to ρ * = 0.5. At densities higher than 0.5, the mixture exhibits the smectic B phase where, as we will analyze later, LCs start to form layers in the direction of the smectic director. Furthermore, in contrast to the 80:20 mixture, we observe that the extent of nematic ordering in the DSS is similar to the LCs in the nematic phase; however, at higher densities, the nematic ordering in the LCs is higher than the DSS [33,34].
We further investigated the relative orientations of the directorsn e andn s for LC and DSS, respectively, in the nematic and smectic phases by calculating the biaxiality parameter, B, defined as follows [34,57]: where the parallel and perpendicular relative orientations of the directorsn e andn s correspond to the value of B equal to 1 and −0.5, respectively. We show the biaxiality parameter with green star symbols in Figure 2. The value of B remains very close to 1, which suggests that the directors of both components of the LC-DSS mixture remain parallel to each other in the nematic and smectic phases. At this point, it is interesting to investigate the degree of magnetization in the LC-DSS mixture. In the 80:20 LC-DSS mixture, where the density of DSS is low, the net magnetization remains negligible in the nematic and smectic phases. However, in the present case, one may expect a ferromagnetic order as the dipolar soft spheres display spontaneous magnetization at high densities [48,58]. In order to investigate the ferromagnetic order in the DSS, we calculated the polar order parameter M = | ∑  These features can be visualized in the snapshots of the mixture shown in Figure 3. The isotropic phase of the LC-DSS mixture, visible at low densities, is characterized by randomly oriented LCs; see Figure 3a. The DSS form chain-like structures that remain randomly oriented and distributed in the system; see Figure 3d. As shown in Figure 3b Figure 3c, the two directors,n e and n s (shown by orange and pink arrows) are aligned parallel to each other.
Afterward, we analyzed the structure of the LC-DSS mixture via the pair correlation functions of different species. As shown in Figure 4, in the isotropic phase, the pair correlation functions for the LCs and DSS show a liquid-like behavior. For both species, the pair correlation function has an oscillatory behavior at a short distance, which decays quickly and the pair correlation functions saturate to 1. The first peak in the g(r) (which accounts for the nearest neighbors) is more pronounced for the DSS than the LCs due to their strong chain-forming tendency.  In the nematic and smectic phases, the projected pair correlation function of both species for distances parallel and perpendicular to the respective directors can be defined as follows [60,61]: δ r ij, /⊥ − r /⊥ θ σ/2 − r ij,⊥/ πNρ * (σ/2) 2 ∆r /⊥ (8) where r ij, =| r ij ·n s,e |, r ij,⊥ =| r ij − r ij ·n e,s n e,s |, and ∆r /⊥ is the volume of the cylindrical shell considered around the particle i. Figure 5a,b shows the projected pair correlation functions for the LCs in the direction parallel and perpendicular ton e , respectively, as a function of distance. In contrast to the isotropic phase, the height of the first peak in both g(r ) and g(r ⊥ ) is more pronounced due to the enhanced density of the system. The location of the first peaks in g(r ) and g(r ⊥ ) are around ∼2.5σ and ∼σ, respectively, indicating the orientational alignment of LCs parallel ton e . However, since the g(r) quickly converges to 1, the layer-wise ordering is not expected.
The longitudinal pair correlation function, g(r ), for the DSS, plotted in Figure 5c, has distinct peaks that occur at an interval of distance σ, confirming the strongly correlated chain-like configuration of the DSS oriented along the directorn s . In the transverse direction, the first peak in the pair correlation function occurs around 2σ, suggesting that the second DSS chain is two-particle diameter apart; see Figure 5d. This further suggests that an LC layer separates two alternating DSS chains, or, as shown in the inset of Figure 5d, the DSS chains remain localized in the cylindrical channels formed by the LCs [36]. However, since the mixture is in the liquid state, the cylindrical channels formed by the LCs are transient, and no long-range order is expected, as revealed by the quick decay of the transverse pair correlation functions for the DSS (see Figure 5d).
In the smectic phase, as shown in Figure 6a, the longitudinal pair correlation function for the LCs has an oscillatory behavior. The peaks in g(r ) appear at a regular interval ∼2.8σ suggesting a layer-wise ordering of LCs along the director,n e . The transverse pair correlation function, shown in Figure 6b, suggests a hexagonal ordering of LCs inside a layer, which would then correspond to the smectic B phase at high densities [33,34,61]. However, the smectic ordering requires further investigation. The situation for the DSS, shown in Figure 6c,d, is very similar to the nematic case. The component of the pair correlation function parallel to the directorn s has distinct peaks at distances comparable to single-particle diameter. The perpendicular component also shows weak oscillations, indicating weak ordering in the transverse directions.

Discussion
To summarize, using extensive molecular dynamics simulations, we have presented a study on the self-assembly of a 50:50 binary mixture of liquid crystals and magnetic nanoparticles. Our results demonstrate that upon compression, the LCs show isotropic, nematic, and smectic B phases at low, intermediate, and high densities, respectively. We identified different phases using nematic order parameters for LCs. The DSS form chains, which remain randomly oriented in the isotropic phase. In the nematic and smectic phases, the DSS chains follow the LCs and orient along with the director of LCs. The orientational ordering of the two species was quantified using nematic order parameters. We found that, in contrast to the 80:20 mixture, the extent of nematic ordering in the DSS in the present case is very similar to the LCs in the nematic phase. However, at high densities, the DSS display lower nematic ordering than the LCs, similar to the 80:20 mixture. The relative alignment of the directors of both components was investigated via the biaxiality parameter, which reveals that in the nematic and smectic phases, the DSS chains follow the average orientation of the LCs. Although the DSS form ferromagnetic chains in the nematic and smectic phases, the macroscopic magnetization of the LC-DSS mixture remains very low, as reflected by the polar order parameter. The mixture does not show demixing or microphase separation in the range of considered densities.
The structure of the LC-DSS mixture was explored using pair correlation functions. In the nematic and smectic phases, the pair correlation functions were resolved in the parallel and perpendicular directions of the respective directors. Our results suggest that albeit with an equal number density, the DSS follow the phases of LCs and do not show ferromagnetic phase at high densities of the LC-DSS mixture. The transverse pair correlation function of the DSS indicates that the two DSS chains remain separated by a column of the LCs. In this context, it will be interesting to explore mixtures with higher number densities of DSS, which may show spontaneous ferromagnetic transitions at high densities. Furthermore, the structure and rheology of the LC-DSS mixture in the presence of an external magnetic field should be interesting as several intriguing effects are expected due to the interplay of the microstructure, external flow, and magnetic fields [51,62].

Data Availability Statement:
The supporting data for this work are available upon request.