Origin and structure of liquid crystalline Blue Phase III

We report here an off-lattice NVT molecular dynamics simulation study of a system of polar chiral ellipsoidal molecules, which spontaneously exhibits Blue Phase III (BPIII), considering coarse-grained attractive-repulsive pair interaction appropriate for anisotropic liquid crystal mesogens. We have observed that suitable selection of chiral and dipolar strengths not only gives rise to thermodynamically stable BPIII but novel Smectic and Bilayered BPIII as well. Further, we have demonstrated that the occurrence of BPIII and its layered counterparts depend crucially on molecular elongation.

origin and structure of liquid crystalline Blue phase iii tanay paul & Jayashree Saha * We report here an off-lattice NVT molecular dynamics simulation study of a system of polar chiral ellipsoidal molecules, which spontaneously exhibits Blue Phase III (BPIII), considering coarsegrained attractive-repulsive pair interaction appropriate for anisotropic liquid crystal mesogens. We have observed that suitable selection of chiral and dipolar strengths not only gives rise to thermodynamically stable BPIII but novel Smectic and Bilayered BPIII as well. Further, we have demonstrated that the occurrence of BPIII and its layered counterparts depend crucially on molecular elongation.
Chiral liquid crystals, in addition to broken rotational and translational (optional) symmetries, belong to the symmetry group that lacks reflection symmetry. Chiral liquid crystal molecules exhibit translationally disordered cholesteric phase where the nematic director precesses about the helical axis. The isotropic-cholesteric phase transition can go through a cascade of intermediate Blue Phases (BP), which have drawn widespread interest in recent technological and biological fields of research. Experimental evidences [1][2][3][4] indicate that though the BPI and the BPII possess cubic structure for the double-twist cylinders, BPIII, the so-called 'blue fog' 3,5 , seems to be amorphous. Apparently, the BPIII does not exhibit Bragg scattering 6 whereas it is attributed with strong optical activity along with thermodynamic stability but is lacking birefringence and having mechanical properties more like isotropic phase. Among chiral phases, the origin and structure of the BPIII phase still remain as a long-standing puzzle, despite extensive endeavour 1 , because experimental explorations have been very much hindered by the occurrence of the BPIII phase in an extremely narrow temperature (or other factors affecting phase transitions) interval and optical effects are manifested only at very short wavelength range.
There are theoretical studies on the structure of the BPIII proposing that it is a quasicrystal 6,7 . Others conclude that BPIII is an amorphous phase having 'spaghetti' like arrangement of double-twist cylinders (DTCs) 8 , or could be a metastable state 3 . Electron micrograph experiments support an amorphous structure and thus eliminate the quasicrystal theory. A recent computer simulation study provides evidence that the structure of BPIII is basically an amorphous network of disclinations 9 . Theoretical 1,7-9 models proposing amorphous 'spaghetti' like arrangement of double-twist cylinders and amorphous network of disclinations helped in determining phase properties. However, a complete understanding of these behaviours requires the information depicting molecular arrangement at the microscopic level responsible for giving rise to the BPIII, which is yet to be achieved. In the present work, we have considered systems composed of molecules interacting through chiral, dipolar and attractive-repulsive van-der-Waals' type interactions. In this Molecular Dynamics simulation study, a microscopic structure resembling BPIII phase and its layered and novel bilayered counterparts have been realized. We also have found that molecular elongation has supported more efficient self-assembly, thus has acted as a stimulating molecular feature to widen stability region. Additionally, the study of the chain-length size 10 of the chiral molecules is important for many technological devices [11][12][13] , lipid bilayer cell-membranes 14  www.nature.com/scientificreports/ The first, second and third terms in Eq. (1) represent the chiral, dipolar and Gay-Berne 21 interactions respectively. Here, r ij is the separation between the centers of mass of the molecules i and j; û i and û j represent the orientations of the molecular long axes of respective molecules with respect to the laboratory axes (Fig. 1b).
Here c is the chirality strength parameter and the handedness of the twist depends on its sign. The form of the chiral interaction potential U C (� r ij ,û i ,û j ) (first term in Eq. (1)) used in our study can directly be obtained from the multipole expansion of electrostatic interaction 22 . To incorporate the ellipsoidal shape of the molecules the orientation dependent well-depth term ǫ(r ij ,û i ,û j ) and separation term R ij = [r ij − σ (r ij ,û i ,û j ) + σ 0 ]/σ 0 have been taken as of Gay-Berne type 21 , where σ (r ij ,û i ,û j ) is the contact distance i.e. the minimum approachable distance between two molecules. σ 0 corresponds to the contact distance at side-by-side ( û i �û j and û i , û j ⊥r ij ) configuration which is also the breadth of the molecule. The term {(û i ×û j ) ·r ij }(û i ·û j ) induces a twist giving rise to chiral phases.
The second term in Eq. (1) is the dipolar interaction term. The polar part of a molecule has been represented by a single terminal point dipole (Fig. 1a), positioned at 0.5σ 0 distance from one terminal point of a model molecule. The orientation of the dipole has been fixed at an angle of 90 • relative to the molecular long axis 19 . Here, � r d ij = r d ijr d ij is the separation between the point dipoles fixed at two molecules i and j (Fig. 1b), û d i and û d j are their respective orientations relative to the simulation box; 'd' suffix corresponds to the dipole. The dipole moment vectors of respective point dipoles are � is the magnitude of the dipole moment in reduced unit, ε 0 is the well-depth for a pair of molecules in side-by-side configuration. The long-range correction of the dipolar interaction has been taken care of by standard Reaction Field technique [23][24][25] . In this technique, a sphere of a cut-off radius r RF is considered around a particular molecular dipole and all the dipoles outside this cut-off sphere are considered to form a dielectric continuum of dielectric constant ǫ RF , which produces a reaction field inside the sphere. The intensity of the reaction field applied on i-th molecule is given by 17 , In our study, the values of cut-off radius r RF and continuous dielectric constant ǫ RF have been taken as r RF = 0.5× the cubic simulation box side-length and ǫ RF = 1.5 26 respectively.
The GB parameters κ [length to breadth ratio], κ ′ [well-depth ratio = ε e /ε 0 ; ε 0 and ε e being the well-depths in the side-by-side and end-to-end ( û i �û j ; û i , û j �r ij ) configurations respectively] and the relative well depth controlling parameters µ and ν have been taken as 3.0, 1/5, 1 and 2 respectively 27 . Potential cut-off radius ( r c ) has been taken as equal to the half of the simulation box-length. To minimize computation time we have considered www.nature.com/scientificreports/ cut-off distances the same for both GB and chiral interactions. However, both the potential functions have been smoothed out at the cut-off boundary by shifting them by the amount U cut = U(r c ) 17 .
We have used a Leap-Frog algorithm 28 for Damped Force method (Hoover's thermostat 29 ) to solve the equations of motions in this NVT-Molecular Dynamics (MD) simulation study 17 which is applicable for the molecules having both translational and rotational degrees of freedom. Scaled density ρ * ρ * ≡ Nσ 3 0 V has been set to 0.30 for κ = 3.0 26,27 . N is the total number of molecules and V is the simulation box volume. For each system, a well equilibrated isotropic phase has been used as the initial configuration and then the scaled temperature ( T * ≡ k B T/ε 0 , k B = Boltzmann constant) has been decreased gradually to study the phase change. For a particular system at each temperature step, the initial configuration has been a previous higher temperature equilibrium phase and simulation run of 10 6 steps has been performed to obtain the equilibrium configuration at that temperature. To meet equilibrium criteria, at each MD step, the average energy of the system has been calculated keeping its variation with MD steps within 2% Root Mean Squared fluctuation about a mean value at equilibrium. Average values have been calculated from the next 10 5 steps after equilibration. As the twisted cylindrical domains remain oriented at random in the BPIII, so they have no biasing in twist directionality, we have used conventional cubic periodic boundary condition. In our MD simulation, the force and torque equations for all the N molecules have been solved through the Leap-frog Verlet integrator method to follow the real trajectory of the system in phase space. For this reason, MD simulation is much expensive computationally 30,31 . To understand the phase structure of polar chiral liquid crystal molecules, we have used an NVT algorithm in this work which is developed by us and used in our earlier works 19,32 .

Results
For characterization of the phases, several distribution functions have been computed which have supported the formation of Blue Phase III (BPIII). To check system size effect we have simulated system sizes corresponding to N = 500, 864, 1372 and some results for N = 2048 . The results have shown qualitatively similar phase sequence behaviour, though stabilization of a particular phase occurs at a different range of scaled temperature, usually shifted to a lower value for larger system size, which has been due to the finite size effects. In this paper we have presented the results of the system sizes N = 1372 and 2048. For all the system sizes with κ = 3.0 , phase properties have been studied for some selected discrete values of the chiral strength parameter c. The value of the reduced dipole moment µ * has been set fixed at a typical value of 1.0 19 . For c = 3.0 and 4.0, Blue phases have been generated from a higher temperature isotropic phase (at T * = 6.0 ) by decreasing temperature. To find out the orientational correlations between the molecules, longitudinal orientational correlation functions 33 S 220 (r * � /d) and S 221 (r * � /d) have been calculated as functions of the intermolecular separation r * � (in units of σ 0 ) measured along a reference axis and further scaled by a distance d which was related to the periodicity of the phase studied www.nature.com/scientificreports/ in Ref. 20 . Here, these functions have been calculated considering minimum image convention and taking the reference axis along one of the simulation box axes and the simulation box-length as the scaling length d. The mathematical forms of these functions are given by, where ... indicates the average over all molecular pairs separated by a distance r * � /d along the chosen reference axis. Here, the function S 220 has a maximum for two parallel molecules, whereas, the function S 221 has an extremum for two side-by-side molecules with 45 • angle between their long axes and thus both help characterizing the chiral phases. The plots of S 220 (r * � /d) and S 221 (r * � /d) (Fig. 2) show qualitatively the same variation as that of a Blue phase, i.e. depending on pitch length they vary approximately as sinusoidal functions with distance 19,20 but the plots are not smooth, because in BPIII cylinders formed from twisted molecular arrangements are themselves twisted and are having comparatively shorter lengths, but in other Blue phases twisted molecular organization makes double-twist cylinders with straight-line symmetry axes spanning over the whole system. www.nature.com/scientificreports/ For other chiral phases obtained with lower c values, these plots were quite smooth and showed sinusoidal type variation 19,20 indicating the presence of definite orientational correlation between the molecules spreading over the whole simulation box. In this case, with relatively higher values of c, the presence of a typical variation of the orientational correlation functions with molecular separation indicates that the obtained phases are not isotropic or nematic, but a chiral one, as isotropic plots show no orientational correlation at all whereas a perfect nematic phase gives almost a straight horizontal line of high correlation value with almost no variation showing a strong parallel orientational correlation between all the molecules. The plots indicate these phases are not like cholesteric, BPI or BPII 19,20 , but have a different Blue phase character originated from chiral molecular organization packed within short size twisted cylindrical arrangements. Snapshots of the configurations of these phases with different c and κ values are presented in Fig. 3. Snapshot of a typical configuration obtained for κ = 3.0 , c = 3.0 is shown in Fig. 3a. Here in the case of BPIII these cylinders themselves are seen twisted in a random fashion. Careful observation of the configurations obtained in these cases reveals that the double-twist cylinders are not extended straight from a box face to the opposite face which occurs in case of other BPs. Snapshots depicting the configurations within five consecutive layers of same thickness perpendicular to one box axis are provided in Fig. 4a-e, where circles are drawn to show the cross-sections of some of the double-twist cylinders. In Fig. 4a-   www.nature.com/scientificreports/ in each case are not straight, but they are tangled, as seen from the configurations. Also, at a higher value of κ , stable equilibrium Blue phase has been realized at a higher value of the scaled temperature T * . As for example, for systems with N = 1372 and c = 3.0 , the typical values of T * at which equilibrium Blue phases have been obtained are T * = 1.8 for κ = 3.0 , T * = 7.5 for κ = 4.0 and T * = 11.0 for κ = 5.0 . For further investigation, the simulation box has been divided into some planar layers perpendicular to one of the simulation box axes and then S 220 (r * � /d) and S 221 (r * � /d) have been calculated for all the planes separately, considering r * � along one of the simulation box axes parallel to the planes. Plots of these functions (Fig. 5) show multiple peaks, not coinciding at same points for all of the layers, but they are at different points for different layers, indicating the fact that the cross-sections of the double-twist cylinders in all planar layers are not at same positions. The number of peaks is more for higher values of c showing lower pitch values (Fig. 6). This indicates that the number of double-twist  www.nature.com/scientificreports/ cylinders increases with the increase in the chiral strength parameter c. These double-twist cylinders are found to be intertwined to form spaghetti-like structures 8 as speculated for BPIII. With further decrease in temperature, the formation of smectic layers has started to develop in addition to the orientational arrangement of the BPIII phase. For a system with a higher value of κ , smectic layers have started forming more efficiently at a relatively higher value of the scaled temperature. As for example, for systems with N = 1372 and c = 3.0 , the typical values of T * at which equilibrium phases with smectic domains have been obtained are T * = 1.4 for κ = 3.0 , T * = 6.5 for κ = 4.0 and T * = 9.5 for κ = 5.0 . Smectic domains are more prominent with higher κ values, i.e. for the systems with higher molecular lengths like κ = 4.0 and 5.0, rather than κ = 3.0 . Snapshots of the configuration obtained with κ = 5.0 are presented in Fig. 7 dividing the simulation box into five planar layers perpendicular to x-axis for visualization of the twist of the double-twist cylinders across the planar layers and formation of layered domains. The plots of the pair distribution function g(r * ) = V � i j� =i δ(r * − r * ij )�/N 2 ( r * is the scaled intermolecular separation = r/σ 0 ) for these phases (Fig. 8), at respective scaled temperature T * at which these smectic phases are stable, show short range positional order without any long range positional order.
Additionally, we have checked the effect of dipole strength on the layer formation 19,32 . The presence of dipoles gives more stability to the smectic layers. Interestingly, for a higher value of µ * novel bilayered smectic BPIII (Fig. 3d, g) has formed, by decreasing temperature from a relatively higher temperature phase. For a system with N = 1372 and c = 3.0 , by decreasing the scaled temperature from T * = 6.5 to T * = 6.0 bilayered arrangement in the smectic domains has formed where κ = 4.0 and µ * = 1.6 (Fig. 3d), whereas in a system with κ = 5.0 and µ * = 1.4 (Fig. 3g) decreasing the temperature from T * = 9.5 to T * = 9.0 bilayered arrangement has formed in the smectic domains. For κ = 3.0 formation of the bilayer is not so prominent. Plots of g(z * ) and g d (z * ) (where the former is the pair correlation function for the molecular centers of mass as a function of the projection ( z * = z/σ 0 ) of separation vector along the axis of a cylindrical domain considered around each molecule and the latter is the same for dipolar positions) have been drawn (Fig. 9)  www.nature.com/scientificreports/ plot. Thus, g(z * ) = �δ(z * − z * ij )�/πR 2 ρ * (where R is the radius of the cylindrical sampling region) and g d (z * ) has been calculated in the same way for the dipolar positions considering cylindrical domain of the same size. The peaks of the function g(z * ) indicate the presence of the molecular centers of mass in planar layers perpendicular to the symmetry axes of the cylindrical sampling regions and similarly, the peaks of the function g d (z * ) indicate the presence of the dipolar positions in similar planar layers. For κ = 4.0 (Fig. 9a, b) cylindrical domain of length 5σ 0 and radius 4σ 0 has been taken. Those values for κ = 5.0 (Fig. 9c, d) are 7σ 0 and 5σ 0 respectively. For both the κ values, the peaks of comparable heights for both the functions occur at nearly same positions for lower values of µ * (Fig. 9a, c) indicating non-bilayer arrangements, but when the value of µ * increases the peaks of comparable heights for g d (z * ) occur alternately with the peaks of g(z * ) (Fig. 9b, d) indicating the presence of small bilayered smectic domains in the BPIII. Due to the smaller size of the smectic domains, we have obtained fewer peaks in these plots. A snapshot of the bilayered BPIII is presented in Fig. 10 by periodic repetition of the simulation box, where bilayered domains can be identified with the dipolar ends of the molecules (shown in black dots) of two adjacent smectic layers clustered together. conclusion Chirality, which is a very important phenomenon occurring in various fields of nature, when added to liquid crystal systems generates many phases featuring fantastic and remarkable properties and having wide applicability in technology, medical science, agriculture and chemical industry. Blue phases are one of this kind. While other liquid crystalline phases have rotational and translational symmetries, chiral Blue Phases (BP), the so-called 'crystalline liquids' , have additional symmetries of conventional solid crystals i.e. cubic crystal symmetries. Being a member of liquid crystal family, this property of BP makes them interesting both as a curious subject for scientific studies and a material of immense potential for industrial use. Theoretical and computer simulation studies on BPI and BPII have proposed that for relatively high chiral strength, the simple helical structure of the cholesteric phase is energetically less stable locally than a 'double-twist cylinder' structure. In a double-twist cylinder, local directors rotate simultaneously about any radius of the cylinder and the local director is parallel to the cylinder symmetry axis at its center. Such double-twist cylinders do not fit in three-dimension to fill the whole space and thus disclinations or defects are formed. Among three types of BP observed experimentally without electric field, in BPI, the arrangement of these disclinations is body center cubic and in BPII is simple cubic. However, BPIII, the so-called 'blue fog' 5 , is amorphous, the structure of which is not very clear till date. Our simulation work supports the theoretical proposition of the 'spaghetti' like arrangement of the double-twist cylinders as a model of BPIII and at the same time provide microscopic description of the molecular arrangement. In this simulation study, we have found that higher chiral strength in the system induces inter-twinning of double-twist cylinders which eventually gives rise to a structure resembling BPIII phase. Additionally, we have found that higher values of dipole moment can induce novel bilayer arrangement in the smectic BPIII. Change in phase properties with www.nature.com/scientificreports/ increase in molecular length has also been studied. Our observation is that the Isotropic-BPIII transition temperature is higher for a system with higher molecular length. Greater the molecular length, the required value of c is lower at which BPIII is formed. We also have observed that molecular elongation favours the efficacious formation of the smectic BPIII, whereas, increased dipolar strength is the key to give rise to novel bilayered BPIII. The focus of our work is to find out the proper contribution of various physical microscopic interactions responsible for the realization of BPIII and its novel layered counterparts. Consequently, this work indicates a way that may help to minimize present experimental hurdle due to significantly narrow stability region 1 of the BPIII phase. We hope, the present coarse-grained simulation study will help in the fundamental qualitative understanding of the molecular level arrangement in the BPIII and the coveted structure-property relationship.