Influence of Volume Fluctuation on End-over-End Rotational Dynamics for a Nematic Liquid-Crystal Phase Molecular-Dynamics Simulation
Molecular simulations by single anisotropic particles have been used to analyse structures and dynamics of liquid crystal systems . Some mesogenic model systems have not only exhibited realistic behaviour of liquid crystals, but have also predicted interesting features or new liquid crystal phases. Further, simulation studies using simple mesogenic models are able to derive data ranging from the molecular to the meso-scale level  because of efficiency in computation, although characteristic details such as conformation and electrostatic properties are inevitably coarse-grained.
The rotational dynamics of liquid crystals in nematic and isotropic phases have been investigated in experiments  and theoretical studies . In particular, end-over-end rotational motion, which can be measured by dielectric relaxation experiments, is one of the characteristic dynamic properties of liquid crystals. The relaxation behaviour for end-over-end rotational motion of liquid crystalline molecules and for the realignment process by an external field in the nematic phase is very slow-mode, due to the orientational order in the long range. Therefore, a long time is required to find the average dynamic property in a molecular simulation.
The concept of retardation factor g|| was introduced by Maier and Saupe [4 - 6] to interpret the relaxation behaviour of nematic-phase molecular rotation. The retardation factor is defined as g||=t||/t0, where t|| and t0 are relaxation times for the end-over-end rotational motion of molecules along the long axis when potential mean torque acts and vanishes, respectively. Recently, Bates and Luckhurst have shown that molecular-dynamics simulation for a mesogenic model system in which the interactions between the molecules are represented by the Gay-Berne potential  GB(4.4,20.0,1,1) at constant volume condition is able to directly estimate the Maier-Saupe strength parameter, which can be measured from dielectric relaxation spectroscopy . They found a density dependence of the retardation factor, in which a free-volume effect may play a role. Further, the relaxation time for end-over-end rotational motion in the nematic phases was found to be larger than the values extrapolated for the isotropic phase at the same temperature. They also mentioned that the Martin-Meier-Saupe solution to the theory underestimated the strength parameter of the Maier-Saupe form [4 - 6] s(@q/RT) by roughly 15%, and the result from the constant volume simulation underestimated the Martin-Meier-Saupe solution . In particular, the discrepancy between the Martin-Meier-Saupe solution and the results in the simulation were significant for highly ordered states , although the generic mesogenic model has shown that the structure and other dynamics behave very realistically in liquid crystalline phases . Our simple question is whether volume fluctuation affects the slow dynamics of rotational motion, for which such a long simulation is required to determine average properties. In fact, the fundamental significance of pressure control in molecular simulation for the liquid crystal systems has been pointed out [11, 12]. The estimation for the slow-process dynamics might then be influenced by the pressure control in a simulation. Since most experiments are performed under a constant pressure condition, knowing the effect of volume fluctuation on the dynamics in a molecular simulation is desired to clarify the method's validity for evaluating real liquid crystal properties. However, there is no report of the influence of volume fluctuation on the end-over-end rotational of molecules in the nematic phase. In this study, we discuss the influence of volume fluctuation in the simulation of nematic-phase rotational motion by a constant pressure molecular-dynamics simulation.
2 Computational Details
A molecular-dynamics simulation at isothermal-isobaric conditions (NPT) corresponding to those used in prior experiments was performed in order to investigate the dynamics of the molecules. Temperature and pressure were controlled using the method of Melchionna et al. , in which the equations are written to generate the correct trajectory for the NPT ensemble. This method is based on the extended system method by Andersen , Nose  and Hoover . As Dominguez et al.  also found, the three components of pressure tensors in the nematic and isotropic phases were the same, although the pressure tensor in the smectic phase was significantly different. A simulation box of cubic shape was used to investigate the dynamics in the nematic and isotropic phases. To solve the equations for the case of the isothermal-isobaric ensemble, an implicit leap-frog scheme was applied. The generic Gay-Berne model, GB(4.4,20.0,1,1) was used in order to facilitate direct comparison with the previous result . The static properties of the generic Gay-Berne model, GB(4.4,20.0,1,1) have been investigated in detail by Bates and Luckhurst . These parameters are closely related to those proposed by Luckhurst and Simmonds  as reproducing the typical mesogenic core of real molecules. The reduced time step, Dt*(@Dt(e0/ms02)1/2), was taken to be 0.001, where m is particle mass, and s0 and e0 are units for length and energy, respectively. The reduced moment of inertia, I|_*(@I|_/ms02), of the particles was set equal to 1.0. The reduced pressure, P*(@Ps03/e0), was set equal to 2.0, at which the system undergoes a stable transition from the isotropic to crystal phase through a nematic phase as the temperature, T*(@kBT/e0) is reduced. The cut-off distance, rc*(@rc/s0), for the potential energy calculation was set equal to 5.5, and a Verlet neighbour list was used in order to improve the efficiency of the simulation. The second-rank orientational order parameter, <P2>, was determined via the Q-tensor
where uia is a unit vector defining the component, a=x, y or z, of the orientation of the symmetry axis for molecule i, and N is the number of molecules. The second-rank orientational order parameter, <P2>, was defined as the maximum eigenvalue obtained by diagonalizing the Q-tensor; its eigenvector is then identified as the director. System sizes of 2,048, 4,000 and 13,500 were used to check the size dependence of the dynamic behaviour. A relatively large system of 8,788 particles was mainly used in this study.
Molecular reorientation in the nematic phase is expressed by the Wigner rotation matrix [18, 19]:
where W is an angle indicating the molecular orientation with respect to the director. One of the first-rank orientational correlation functions (L=L'=1) for the end-over-end motion of mesogenic molecules with a longitudinal electric dipole is related to the dielectric relaxation experiment . The relevant correlation function is
Figure 1. The first-rank orientational time correlation functions for the system with N of 8,788 in the nematic phase at T* of 1.50, 1.575 and 1.65.
where b is the angle between the molecule and director. In this study, the following formula is used to calculate the correlation function considering the fluctuation of the director,
where n(t) is the director at time t. Configuration data at each of 500.0 time steps for 600.0 t* to 1000.0 t* (0.6 to 1.0×106 time steps) were used to obtain the relaxation decay for end-over-end rotation. Several data sets of different initial configurations were used for analysis of the rotational dynamics.
3 Results and Discussion
The system studied has been known to take isotropic, nematic, smectic A, smectic B and crystal B phases at the reduced pressure P* of 2.0 as the temperature T* is reduced . In this system, the nematic phase was formed in the temperature range of 1.65 to 1.475. The transition behaviour in our simulation was confirmed to be in fair agreement with that of the previous study with the same model system . We also confirmed that the three components of pressure tensor had the same values in the nematic and isotropic phase in the simulations.
The first-rank orientational time correlation functions c||(1)(t) at various temperatures in the nematic phase are shown in Figure 1. The correlation function clearly changed with temperature T* and the second-rank orientational order parameter <P2>, and these functions decayed exponentially at all temperatures in the nematic and isotropic phases. As shown in the next figure, all data could be fitted by single exponential function as follows: c||(1)(t) = c||(1)(0)exp(-t/t||). Here c||(1)(0) is the first-rank orientational correlation function at time 0.
Figure 2. Arrhenius plot of the relaxation time for the end-over-end rotational motion in the isotropic and nematic phases for the system at the reduced pressure P* of 2.0. The open circles indicate the values for the system with N of 2,048. The closed circles, triangles and squares indicate the values for the system with N of 4,000, 8,788 and 13,500, respectively.
Figure 2 shows the Arrhenius plot for the relaxation time t|| of the first-rank rotational correlation function c||(1)(t) in the nematic phase at a P* of 2.00 for the system. As can been seen, the Arrhenius relationship clearly holds. The point at the highest nematic temperature T* of 1.65 (1/T* of 0.606) should be noted. This point slightly deviates from the nematic-phase relation. Further, a few points in the higher region of the isotropic phase also deviate from the straight line. The dynamics for these points are slightly different from those of the other points because of the transition between the nematic and isotropic phases. Here we shall consider only the points at which the Arrhenius relation is valid to discuss the retardation factor g||(@t||/t0). The relaxation time t0 for a hypothetical state in which the potential of mean torque vanishes was obtained by extrapolating the relaxation time for the isotropic phase into the nematic phase assuming the same Arrhenius-type behaviour as in the experiments. As shown in Figure 2, the behaviour for the smaller systems of N = 2,048 was different from the larger systems throughout the temperature range. However, for the larger systems of N > 4,000, there was no obvious difference in the relaxation times in the nematic and isotropic phase. The behaviour for the system with 8,788 was clearly consistent with the values for the larger system of N = 13,500.
In liquid crystalline phases, not only the rotational motion, but also the translational motion is influenced by orientational order. In order to examine the translational dynamics in the nematic phase, the temperature dependence of the calculated translational diffusion coefficients is shown in Figure 3. Here the translational diffusion coefficients of parallel (D||) and perpendicular (D|_) components were evaluated by the following autocorrelation functions:
where v||i(t) and v|_i(t) are parallel and perpendicular components of the linear velocity for molecule i at time t. The data for 5.0t* were used for determining the diffusion coefficients. In this figure, the data for the parallel component of the D|| coefficients at higher temperatures T* of 1.65 and 1.625 (1/T* of 0.606 and 0.615) are not included because they clearly deviate below the line of the Arrhenius relation. The isotropic structure is partially contaminated at these points because of the vicinity of the transition temperature. Similarly, the point at the lowest point T* of 1.475 (1/T* of 0.678) was also not on the straight line. A smectic structure might form locally at the lowest temperature in terms of the dynamics for translational motion although no tendency for the coexistence of smectic and nematic phases is shown in the Arrhenius plot for the first-rank orientational correlation function c||(1)(t) in the nematic phase (Figure 2). Of course, parallel-direction translation motion in the smectic phase is hindered because ofthe layer structure.
Figure 3. Arrhenius plot of the translational diffusion coefficients in the directions parallel to the director, D||, and perpendicular to it, D|_, in the nematic phase for the system at the reduced pressure P* of 2.0 for the system with N of 8,788.
Figure 4 shows the retardation factor g|| in the nematic phase as a function of the Maier-Saupe strength parameter s(@q/RT) obtained from the simulation at constant pressure for the system of N = 8,788 together with two theoretical lines [4 - 6, 9] and the data obtained from the simulation at constant volume . To our knowledge, experimental data on the retardation factor as a function of the second-rank orientational order parameter <P2> are not available, although some data on the retardation factor as a function of temperature are. [21, 22] However, the mean field theory, which was successively described for the nematic phase by Maier and Saupe, can be used to explain the retardation factor g|| as a function of the second-rank orientational order parameter <P2> and/or the mean field-strength parameter, namely the Maier-Saupe strength parameter s. Therefore, we discuss the end-over-end rotational molecular motion in the nematic phase using the order parameter and the Maier-Saupe strength parameter s in this paper. Clearly, the original theoretical curve by Meier and Saupe overestimates g|| compared to the later work , which gives a solution by Martin et al. . Moreover, the results from the simulation at constant volume underestimate g|| as compared with the work by Coffey et al. in the range of all values for s. On the other hand, the values in our simulation at constant pressure show a slight overestimation in the region of s < 6.4 and an underestimation in the region of s > 6.4, which is a highly ordered state at a <P2> greater than 0.75. Therefore, the values from the simulation with pressure control are not in accordance with those from the simulation at constant volume, although the values of the number density r*(@N/V) are very close. More or less, the results from the simulation with pressure control gave a better solution than the results from the simulation without volume fluctuation for the entire range of s studied. The values of number density fluctuated in our simulation because the pressure was kept constant. The values of the density r* changed from 0.170 to 0.183 at a reduction of T* from 1.65 to 1.475. Therefore, volume fluctuation is one possible cause of the disagreement in taking the first-rank rotational correlation function which requires a long time to average. The change in the volume between the maximum and minimum of the simulation box was 1.7% in the nematic phase. This volume fluctuation may not have as much effect in the short term as the second-rank orientational correlation function c||(2)(t). However, this fluctuation might have an effect in estimating a property obtained from a long-term average such as a relaxation time for end-over-end rotational motion. The values of g || at the highly ordered states of <P2> greater than 0.75 show large differences in spite of similar densities.
Therefore, the result from our simulation at constant pressure more closely reproduced the analytical expression by Coffey et al.  throughout a wide range of <P2> and s. In particular, the values for the highly ordered states were improved over the values of the Martin-Meier-Saupe solution by simulation with pressure control. Consequently, the end-over-end rotational dynamics in the nematic phase is sensitive to the ensemble in the molecular-dynamics simulation for the Gay-Berne mesogenic model GB(4.4,20.0,1,1). This means that a fluctuation of volume would be required in the simulation for analysis of slow dynamics in the nematic phase. Finally, in highly ordered states with <P2> greater than 0.75, the values in our study underestimated the Martin-Meier-Saupe solution although they were better than the results in the results previously obtained from a simulation without volume fluctuation . This might be caused by the Gay-Berne model having a simple ellipsoidal shape without flexible parts unlike real mesogens. Thus, considering a realistic complexity such as alkyl chains and polar parts may produce a more concordant result in the highly ordered state.
Figure 4. Retardation factor g|| as a function of the Maier-Saupe strength parameter s(@q/RT) and the second rank orientational order parameter <P2> for a nematic phase. The ordinary components for the relaxation times t0 are estimated from the extrapolation line for the data in the isotropic phase by assuming Arrhenius behaviour. The closed circles indicate the values obtained from the simulation in this study. The open circles show the values obtained from the molecular dynamics simulation at constant volume (r* of 0.17 to 0.19) . The dotted and solid lines are the theoretical curves by Meier-Saupe and Martin-Meier-Saupe solution, respectively.
The influence of volume fluctuation on end-over-end rotational motion in the nematic phase has been discussed by isobaric-isothermal molecular-dynamics simulation for the Gay-Berne mesogenic system. The result obtained from the first-rank orientational correlation function c||(1)(t) was compared with the previous result obtained from a molecular-dynamics simulation for the same Gay-Berne model at constant volume and the Meier-Saupe theory. Arrhenius relations of the relaxation time for the end-over-end rotational motion and of translational diffusion coefficients proved to be significantly valid in the middle of the nematic phase for the system under isobaric conditions. In particular, the volume fluctuation in the simulation improved the estimation for the retardation factor g||. It can thus be concluded that pressure control is essentially required in a molecular-dynamics simulation for the analysis of very slow end-over-end rotational dynamics.
[ 1] P. Pasini and C. Zannoni, Advances in the computer simulations of liquid crystals, Kluwer academic publishers (1998).
[ 2] G. Williams, The Molecular Dynamics of Liquid Crystals, G.R. Luckhurst and C.A. Veracini, Kluwer, Dordrecht (1994).
[ 3] P.L. Nordio, R. Rigatti, U. Segre, Mol. Phys., 25, 129 (1973).
[ 4] W. Maier, A. Saupe, Z. Naturforsch, 13a, 156 (1958).
[ 5] W. Maier, A. Saupe, Z. Naturforsch, 14a, 882 (1959).
[ 6] W. Maier, A. Saupe, Z. Naturforsch, 15a, 287 (1960).
[ 7] J.G. Gay, B.J. Berne, J. Chem. Phys., 74, 3316 (1981).
[ 8] M.A. Bates, G.R. Luckhurst, Mol. Phys., 99, 1365 (2001).
[ 9] A.J. Martin, G. Meier, A. Saupe, Symp. Faraday Soc., 5, 119 (1971).
 M.A. Bates, G.R. Luckhurst, J. Chem. Phys., 110, 7087 (1999).
 H. Dominguez, E. Velasco, J. Alejandre, Mol. Phys., 100, 2739 (2002).
 K.M. Aoki, M. Yoneya, H. Yokoyama, J. Chem. Phys., 118, 9926 (2003).
 S. Melchionna, G. Ciccotti, B.L. Holian, Mol. Phys., 78, 533 (1993).
 H. C. Andersen, J. Chem. Phys., 72, 2384 (1980).
 S. Nose, Mol. Phys., 52, 255 (1984).
 W.G. Hoover, Phys. Rev. A, 31, 1695 (1985).
 G.R. Luckhurst, P.S.J. Simmonds, Mol. Phys., 80, 233 (1993).
 C. Zannoni and M. Guerra, Mol. Phys., 44, 849 (1981).
 C. Zannoni, The Molecular Dynamics of Liquid Crystals, G.R. Luckhurst and C.A. Veracini, Kluwer, Dordrecht (1994).
 W.T. Coffey, D.S. F. Crothers, Yu P. Kalmikov, J.T. Waldoron, Physica A, 213, 551 (1995).
 G. Williams, The Molecular Dynamics of Liquid Crystals, G. R. Luckhurst and C. A. Veracini, Kluwer, Dordrecht (1994), and reference therein.
 R. Righini, Physical Properties of Liquid Crystals: Nematic, D. A. Dunmur, A. Fukuda, and G. R. Luckhurst, Inspec, London (2001), and reference therein.