Molecular dynamics simulation of orientational glass formation in anisotropic particle systems in three dimensions
We propose a simple microscopic model of molecular dynamics simulation to study orientational glass in three dimensions. We present simulation results for mixtures of mildly anisotropic particles and spherical impurities. We realize fcc solids without orientational order in a rotator phase. As the temperature is lowered, the disordered matrix is gradually replaced by four kinds of orientationally ordered, rhombohedral domains. Two-phase coexistence is realized in a temperature window. The impurities serve to anchor the orientations of the surrounding anisotropic particles, resulting in finely divided domains or medium long-range orientational order. We examine the rotational dynamics of the molecular orientations which is slowed down at low . We predict the shape memory effect under a stretching cycle due to inter-variant transformation.
Solid-solid transition Defects and impurities in crystals Glasses
Nonspherical molecules such as KCN can form a crystal without long-range orientational order for mild molecular anisotropy , while liquid crystal phases can appear for large molecular anisotropy. Such crystals are called plastic crystals in a rotator phase. They undergo an orientational phase transition as the temperature is further lowered, where the crystal structure is cubic at high and non-cubic at low . With inclusion of impurities in such solids, the so-called orientational glass has been realized . Around the transitions, a peak in the specific heat  and softening of the shear modulus [1, 2] have been observed. In real systems, the molecules often have dipolar moments, yielding dielectric anomaly. As a similar example, metallic ferroelectric glass, called relaxor, has been studied extensively , where frozen polar nanodomains have been observed.
For one-component anisotropic particle systems, there have been a number of simulations on the statics [4, 5, 6, 7, 8, 9, 10] and dynamics [11, 12] of the orientational phase transition. For two-component anisotropic particle systems, Chong et al. examined the slowing-down of the orientational time-correlation functions around the glass transition . In this paper, setting up a simple microscopic model, we will investigate the formation of orientational glass. In particular, we will examine how impurities can microscopically produce orientational disorder, which has remained unexplored in the literature. Furthermore, alignments of anisotropic particles formig a crystal give rise to lattice deformations. As a unique feature in orientational glass, heterogeneous strains should emerge on mesoscopic scales. In such situations, we may predict soft elasticity and a shape memory effect against applied stress.
We consider a binary mixture of anisotropic particles with number and spherical particles with number . In this paper, we set . The composition of the second species is defined by
We assume relatively small , so the spherical particles may be treated as impurities. The particle positions are written as (). The anisotropic particles are assumed to be axisymmetric; then, their orientation may be expressed in terms of the solid angles and () as The particle sizes of the two species are characterized by lengths and . The pair potential between particles and () is a truncated modified Lennard-Jones potential depending on the particle distance and the directions and . For it is zero, while for it reads
where is the characteristic interaction energy and
The particle anisotropy is taken into account by the angle factor , which depends on the relative direction and the orientations and . In this paper, we assume the following form,
where is the anisotropy strength of repulsion. The () is equal to 1 for ) and 0 for (). In the right hand side, the first (second) term is nonvanishing only when () belongs to the first species. The ensures the continuity of at .
In our system, the total potential energy is the sum , while the total kinetic energy arises from the translational velocity and the rotational velocity as
where and are the masses, and is the moment of inertia of the first species. We set in our simulation. The molecular rotation around the symmetry axis parallel to does not change , so we may neglect its kinetic energy. The Newton equations of motion for translation and rotation are written as 
where and . Hereafter, we measure space, time, and temperature in units of , and , respectively.
Treating equilibrium or at least nearly steady states, we furthermore attach a Nos-Hoover thermostat  to all the particles. That is, we added the thermostat terms and in the right hand sides of Eqs.(6) and (7), respectively, where obeys
with being a short relaxation time taken to be 0.1.
We may envisage the anisotropic particles as spheroids depending on the parameter from Eq.(4). For particles and of the first species, minimization of with respect to yields . For , the short and long diameters, and , are given by
for the perpendicular and parallel orientations of and with respect to , respectively. The inertia momentum of the first species in Eq.(7) is given by
In our simulation, we set and . From Eq.(9) the aspect ratio is . For these mild parameter values, we realized fcc plastic solids without long-range orientational order at relatively high . For large , say 10, we realized liquid crystal phases in our model. It is worth noting that our angle-dependent potential (2) is analogous to the Gay-Berne potential for anisotropic molecules used to simulate mesophases of liquid crystals . Similar angle-dependent potentials have also been used for lipids forming membranes. [17, 18].
3 Four-variant states at fixed volume
Starting with liquid states at , we quenched the system to below the melting and waited for to realize fcc plastic crystals. Here, a few stacking faults were formed in most runs. We then cooled the system to a final temperature and waited for , where orientational order developed. In Figs.1-7, we fixed the system volume and shape imposing the periodic boundary condition. In terms of the molecular volumes and , the packing fraction is given by
We set in Figs.1-7. Then the system length is about and the pressure is about .
In Fig.1, the orientational domains can be seen at , where the particle positions and orientations are nearly frozen in time. For , the system undergoes a structural phase transition, where the lattice structure changes from fcc to rhombohedral almost without dilation change. The lattice constant is equal to through the transition. The ordered phase consists of four rhombohedral variants where the angles of the rhombuses are and . The separation distance of the closely packed planes is increased from to 1.08, because the orientation vectors are aligned in the direction. With increasing , the orientational disorder is gradually increased, leading to pinning of finely divided domains. We eventually obtain a glass state at .
In Fig.2, the anisotropic particles are aligned in the perpendicular directions ) around one or two impurities. This parallel anchoring disturbs the orientation order. In Fig.3, in a (111) plane, we display the four rhombohedral variants with different orientations for . The interfaces between the variants tend to be trapped at clustered impurities. We also notice that the junction angles, at which two or more domain boundaries intersect, are approximately multiples of . Similar patterns were observed on hexagonal planes in a number of experiments on alloys after structural phase transitions .
4 Impurity clustering
Figure 4 shows significant impurity clustering, which appeared during quenching from liquid. The average impurity number per cluster is 5.5 for and there is a big percolated cluster for . Here, two impurities are treated to belong to the same cluster if their distance is smaller than 1.3. In the snapshot of in Fig.1, the impurity clustering is closely related to the heterogeneity in the orientations.
To explain this effect, we compare the solvation energy of a single impurity and that of two associated impurities (see Fig.2). In terms of the potential energy of particle , it is estimated as
Here, is the contribution from the impurities under consideration, the summation is over the nearby anisotropic particles with ( and ), and is the average potential energy of the non-neighbor anisotropic particles separated from any impurities longer than . At , is calculated as for a single impurity and for a dimer, where . The difference is the association energy. At , we have . The total potential energy is thus lowered with impurity clustering in the present case.
5 Coarse-grained orientation order parameter
For each particle of the first species, we introduce the orientation tensor () as
where is the unit tensor. The summation is over the bonded particles of the first species with and is the number of these particles. If a fcc lattice is formed, it includes the nearest neighbor particles, so . We define the amplitude of the orientational order for each anisotropic particle as
Here, in disordered regions due to the thermal fluctuations, but it increases up to unity within rhombohedral domains at low . In Fig.5, the average represents the overall degree of orientation order. In our simulation, we realized only uniaxial states, where we have in terms of the amplitude and the director .
6 Coexistence of high and low temperature phases
In Fig.5, the degree of orientational order increases continuously as is lowered. For small at fixed volume and shape, however, we find coexistence of cubic and rhombohedral regions in a temperature window, which is roughly given by for and for . To examine this coexistence, we introduce the distribution of the orientational amplitude,
Taking the average over six runs we obtained in Fig.6, which exhibits two peaks for at . We also give a snapshot of the ordered anisotropic particles with . For , the interfaces between the ordered and disordered regions are temporally fluctuating with small amplitudes. Small amounts of impurities can pin the interfaces and the ordered regions are stabilized far from the impurity clusters.
7 Rotational dynamics
In crystal and glass, elongated particles often undergo turnover motions, . These events take place in a microscopic time ), while the characteristic time between successive turnover motions grows at low . Let us consider the following angle relaxation function,
where the average is taken over the initial time and over six runs in this paper. Furthermore, in terms of the -th order Legendre polynomials , we may define the -th order rotational relaxation functions as [11, 12, 13]
Here decays on the timescale of , while is unchanged by the turnover motions. In orientational glass, the ultimate relaxation of is due to the orientational configuration changes involving the surrounding particles, so its relaxation time much exceeds .
In Fig.7, we plot and vs for . For they relax considerably due to the thermal motions. The decay of is slower for lower , while apparently tends to a positive constant for large , as in the previous simulation . In our case, this plateau appears because the anchoring of the anisotropic particles around the impurities becomes nearly permanent at low . Thus increases with lowering . In Fig.7, we also present the time-evolution of the relaxation function for , and 20000, which is peaked at . The peak height at increases in time from 0. Also shown is the turnover probability defined by
where we set . We notice that grows linearly in time in the very early stage. Thus we may define from the following linear form,
which holds for . For , the probability of multiple turnovers becomes appreciable resulting in a deviation from the linear behavior (19). However, the scaling form fairly holds for . For , the orientational structural relaxation comes into play and the scaling in terms of does not hold.
8 Shape memory effect
In the presence of multi-variant orientational order, the shape-memory effect (mechanical hysteresis) emerges almost without dislocation formation. This effect is well-known for shape memory alloys such as Ti-Ni , where a structural phase transition is caused by the relative atomic displacements in each unit cell. In our case, under stretching along the axis, domains with nearly parallel to the axis grow yielding an increase in the strain. We are not aware of any experiments on this effect in anisotropic particle systems.
In this mechanical effect, the system shape changes. Thus we assumed a Parrinello-Rahman barostat [14, 21] as well as the Nos-Hoover thermostat. Before streching we prepared a multi-variant initial state at , where the packing fraction in Eq.(11) was 0.75 and the pressure was zero. For we applied the stress setting , where are the stress components with being the space average. We allowed the system to take a rectangular shape. In terms of the system length along the axis, the average strain is
We define effective Young’s modulus by . It follows the effective shear modulus from the formula in classical elasticity, where is the bulk modulus. In our case is about and is much larger than , so . Hereafter, and will be measured in units of .
In Fig.8, we first increased at from 0 up to 0.5 at . In the left panel, we find for , for , and for . In the second range, the solid is soft because the variants elongated along the axis grow due to inter-variant transformation. Here, the orientation vectors of the four variants are roughly given by (a) (0.3,-0.1,0.9), (b) (-0.7,0.1,0.7), (c) (0.7,0.7,0.2), and (d) (0.5,-0.9,0.1) in Fig.8. In the right panel, the variant (a) grows, the variant (b) remains nonvanishing, and the variants (c) and (d) disappear. Secondly, we decreased back to 0 at . In this return path, was about , the disfavored variants no more appeared, and a remnant strain of order remained. Thirdly, at , we increased from 0.05 to 0.2 into the disordered phase. After this heating, the system became isotropic.
We have siudied the orientational glass using the potential (2). Due to the anisotropic factor , the anisotropic particles have the aspect ratio . In our simulation, the anisotropy parameter is 1.2 and the size ratio is 1.1. As a result, fcc plastic crystals appear at relatively high , which undergo the orientation transition at lower . The impurities serve to pin the orientations of the surrounding anisotropic particles.
We briefly summarize our main results. In Fig.1, we have illustrated four-variant ordered states for various at . With increasing , rhombohedral domains are finely divided and their typical size is decreased. In Fig.4, we have shown that the impurities are heterogeneously distributed. In Fig.5, we have plotted the average orientational amplitude , which increases continuously with lowering . In Fig.6, we have shown large-scale coexistence of the disordered and ordered phases in a temperature window. In Fig.7, the orientational time-correlation functions have been plotted for , where decays due to turnover motions. In Fig.8, we have examined the shape memory effect in orientational glass.
We next make some remarks below.
(i) For small , we need to know how the characteristic domain size is determined. For moderate impurity concentrations ( in this work), the orientational disorder is proliferated. These aspects should be further studied.
(ii) The phase transition delicately depends on whether simulations are performed at fixed volume or at fixed stress. Some salient features at fixed stress are as follows. (a) For , a first-order phase transition occurs with a shape change. (b) For , both multi-domain states and single-domain states can be realized for the same parameters as in Fig.8. (c) For , the two-phase coexistence in Fig.6 can still be realized even at fixed stress.
(iii) In our simulation, the system keeps the crystalline order. However, for larger size ratio , say 1.4, we found an increase in the positional disorder leading to polycrystal and positional glass. Competition of positional and orientational glass transitions can then be studied. For large anisotropic parameter , we may also study complicated phase behavior of two-component liquid crystals.
(iv) In this paper, the impurities are spherical and slightly larger than the anisotropic particles. By modifying the potential form, we may also treat other types of impurities. For example, their attractive interaction with the host anisotropic particles may be anisotropic; then, the anchoring can be homeotropic.
Acknowledgements.This work was supported by Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan. K.T. was supported by the Japan Society for Promotion of Science. The numerical calculations were carried out on SR16000 at YITP in Kyoto -University.
-  \NameHöchli U. T., Knorr K. Loidl A. \REVIEWAdv. Phys.391990405.
-  \NameLynden-Bell R. M. Michel K. H. \REVIEWRev, Mod. Phys.661994721.
-  \NameVugmeister B. E. Glinchuk M. D. \REVIEWRev. Mod. Phys.621990993; \NameHirota K., Wakimoto S. Cox D. E. \REVIEWJ. Phys. Soc. Jpn.752006111006.
-  \NameFrenkel D. Mulder B. M. \REVIEWMol. Phys.5519851171; \NameVeerman J. A. C. Frenkel D. \REVIEWPhys. Rev. A4119903237; \NameBolhuis P. Frenkel D. \REVIEWJ. Chem. Phys.1061997666.
-  \NameAllen M. P. Imbierski A. A. \REVIEWMol. Phys.601987453.
-  \NameSinger S. J. Mumaugh R. \REVIEWJ. Chem. Phys.9319901278.
-  \NameVega C., Paras E. P. A. Monson P. A. \REVIEWJ. Chem. Phys.9619929060; \SAME9719928543.
-  \NameRadu M., Pfleiderer P. Schilling T. \REVIEWJ. Chem. Phys.1312009164513.
-  \NameMcGrother S. C., Williamson D. C. Jackson G. \REVIEWJ. Chem. Phys.10419966755.
-  \NameMarechal M. Dijkstra M. \REVIEWPhys. Rev. E772008061405.
-  \NameDe Michele C., Schilling R. Sciortino F. \REVIEWPhys. Rev. Lett.982007265702.
-  \NameCaballero N. B., Zuriaga M., Carignano M. Serra P. \REVIEWJ. Chem. Phys.1362012094515.
-  \NameChong S.-H., Moreno A. J., Sciortino F. Kob W. \REVIEWPhys. Rev. Lett.942005215701; \NameChong S.-H. Kob W. \REVIEWPhys. Rev. Lett.1022009025702.
-  \NameAllen M. P. Tildesley D. J. \BookComputer Simulation of Liquids \PublClarendon Press, Oxford \Year1987.
-  \NameNosé S. \REVIEWMol. Phys.521984255.
-  \NameGay J. G. Berne B. J. \REVIEWJ. Chem. Phys.7419813316; \NameBrown J. T., Allen M. P., del Rio E. M. de Miguel E. \REVIEWPhys. Rev. E5719986685.
-  \NameDrouffe J.-M., Maggs A. C. Leibler S. \REVIEWScience25419911353.
-  \NameNoguchi H. \REVIEWJ. Chem. Phys.1342011055101.
-  \NameKitano Y., Kifune K. Komura Y. \REVIEWJ. Phys. (Paris)491988C5-201; \NameManolikas C. Amelinckx S. \REVIEWPhys. Stat. Sol. (a)601980607.
-  \NameSarkar S., Ren X. Otsuka K. \REVIEWPhys. Rev. Lett.952005205702; \NameWang Y., Ren X. Otsuka K. \REVIEWPhys. Rev. Lett.972006225703.
-  \NameParrinello M. Rahman A. \REVIEWJ. Appl. Phys.5219817182.