Dynamic phenomena and viscous properties in a liquid crystal: A theoretical treatment and molecular dynamic simulations

Abstract

The molecular dynamics (MD) simulations, based on a realistic atom–atom interaction potentials were performed on 4-heptyloxy-4′-cyanobiphenyl (7OCB) and 4-hexyloxy-benzylidene-4′-amino-benzonitrile (HBAB) in nematic phase. The set of the order parameters (OPs) S_2L (L = 1, 2, 3), rotational self-diffusion (RSD) coefficient D_⊥, rotational γ_i (i = 1, 2) and Leslie α_i (i = 1, … , 6) viscosity coefficients, the set of the orientational correlation times ${\tau }_{j0}^i$ (i = 1, 2; j = 0, 1) for 7OCB and HBAB in the nematic phase are calculated. Reasonable agreement between the calculated and the experimentally obtained data on S₂, ${\tau }_{00}^1$, and γ₁ has been obtained.

Introduction

A key research opportunity for liquid crystals (LCs) is their potential use in LC displays and nonlinear optical ultra-fast integrated electronic electrooptic modulators needed for ultra-fast (>10 GHz) optical switching in telecommunications and computing [1]. During the last decade, there has been a significant amount of research undertaken with the aim to decrease the switching time of the LC devices [2]. This means that the rotational viscosity and self-diffusion phenomena play a crucial role in the processes responsible for faster switching. But the theoretical treatment of dynamical processes of flexible molecules in an anisotropic medium is not an easy task [3], [4]. In contrast to isotropic liquids, where a single rotational viscosity coefficient (RVC) is sufficient to characterize the flow properties, LCs and other anisotropic systems require several such coefficients. A common approach for determination of viscosity coefficients uses statistical–mechanical theory [5], [6], [7], [8] and is based on the rotational diffusion model [3]. This model, in turn, assumes that the reorientation of the molecules in LCs can be described by a stochastic Brownian process for molecular reorientations in which each molecule moves in time as a sequence of small angular steps caused by collisions with its surrounding molecules and under the influence of a potential of mean torque set up by these molecules. Each molecule is characterized by a rotational diffusion tensor whose principal elements (D_xx = D_yy = D_⊥, D_zz = D_‖) are determined in a frame fixed on the molecule. In such a treatment, the system is determined by the time-dependent single-particle orientational distribution function (ODF) governed by a Fokker–Planck-type kinetic equation. The validity of various theories where the RVCs were expressed as certain polynomials of equilibrium orientational order parameters (OPs) and inversely proportional to the rotational diffusion coefficient (RDC) in the nematic phase [5], [6], [7], [8] has been carried out [9], [10], [11], [12], [13], [14], [15].

In the present study, we extend of our previous investigation of the viscous properties of 4-pentyl-4-cyanobiphenyl (5CB) [9], [10], [15] to two new compounds, 4-heptyloxy-4′-cyanobiphenyl (7OCB) and 4-hexyloxy-benzylidene-4′-amino-benzonitrile (HBAB). The RVCs will be determined using the orientational time correlation functions (TCFs). These functions will be obtained from the MD trajectory.

Measurements of the viscous coefficients of nematic LCs (NLCs) have been a subject of interest from the beginning of the investigations of these materials [16], [17], and continue to do so. The non-equilibrium molecular dynamics (MD) simulations have been used to study the Miesowicz and Leslie viscosities as functions of density and temperature [18], [19], [20]. Using the Gay-Berne potential Sarman [19] and Wu et al. [20] also performed non-equilibrium MD simulations of LC shear flow and calculated the shear and rotational viscosities in the nematic phases.

The reason why many reports are still concerned with the dynamical processes of flexible molecules in an anisotropic medium is that the measurements of motional constants are far from being trivial and these constants play a crucial role in the processes responsible for fast switching of the LC devices.

The classical approach to the viscosity of LCs, in the framework of the Ericksen–Leslie theory [21], [22], assumes that the flow properties can be described by six Leslie coefficients α_i (i = 1, … , 6) which satisfy the Parodi relation α₂ + α₃ = α₆ − α₅. Hence only five of the six coefficients α_i are independent. In the isotropic phase all α_i, except α₄, equal zero, and α₄ is the shear viscosity of an isotropic liquid. In the LC phase the coefficient γ₁ and γ₂ are called the rotational viscosity coefficients. They determine, in absence of hydrodynamics flow, the dissipation of energy $\frac{1}{2}{\gamma }_1{\left(\frac{d\hat{\mathbf{n}}}{\mathit{dt}}\right)}^2$ due to only the rotation of the director $\hat{\mathbf{n}}$ under the action of external forces.

A number of statistical–mechanical approaches (SMAs) for theoretical treatment of rotational viscosity [5], [6], [7], [8], [9], [10], [15] have been proposed. These theories are based on the rotational diffusion model [23]. These approaches rest also on the concept of treating of the phenomenological stress tensor (ST) $\overline{\sigma }$ as an average of its microscopic equivalent σ based on an appropriate nonequilibrium ODF f(t, u) in the form $\overline{\sigma }=\int d\mathbf{u}f\left(t\text{,}\mathbf{u}\right)\sigma$, where u is the unit vector along the molecular symmetry axis. Under the assumption that microscopic molecular motion in the NLC can be considered as a rotational Brownian motion in an external mean potential U, the following kinetic Fokker–Planck equation determines the ODF [10], [15]

$$\frac{\mathit{df}\left(t\text{,}\mathbf{u}\right)}{\mathit{dt}}-{\nabla }_{\mathbf{u}}·\left[{\mathit{fD}}_{\perp }{\nabla }_{\mathbf{u}}\left(\lg f+\frac{U}{k_{B}T}\right)\right]+{\nabla }_{\mathbf{u}}·\left[f\mathbf{M}·\mathbf{u}+\mathit{pf}\left(\mathcal{E}-\hat{\mathbf{n}}\hat{\mathbf{n}}\right)·\mathbf{W}·\mathbf{u}\right]=0\text{,}$$

where k_B is the Boltzmann constant, T is the temperature, p =  (a² − 1)/(a² + 1), a is the length-to-breadth ratio of the molecule, $\mathcal{E}$ is an unit tensor, and M and W are the symmetric and antisymmetric parts of the flow velocity gradient [21], [22] $v_{i\text{,}j}=\frac{\partial v_j}{\partial x_i}$. The external mean potential in Eq. (1) given by the form $\frac{U\left({\theta }_i\right)}{k_{\mathrm{B}}T}=S_2{P}_2\left(\cos {\theta }_i\right)=S_2\left(\frac{3}{2}{\cos }^2{\theta }_i-\frac{1}{2}\right)$, where θ_i is the polar angle between the director and the molecular symmetry axis, $S_{2L}=\int P_2\left(\cos {\theta }_i\right)f_0\left(\cos {\theta }_i\right)d\left(\cos {\theta }_i\right)$ is the even-rank OPs, P₂(cos θ_i) denotes the Legendre polynomials, and f₀(cosθ_i) is the local equilibrium ODF. Here, we employ a small gradient approximation for f = f₀(1 + h), where the stationary correction ∣h∣ ≪ 1.

The average of the symmetric part of ST ${\overline{\sigma }}^s$ can be calculated exactly without solving the Fokker–Planck equation, and by just replacing the statistical–mechanical averaging with nonequilibrium distribution function by the expression containing only the known equilibrium one [5], [15]. As a result, a set of microscopic expressions for the symmetric Leslie coefficients and γ₂ can be written as [5]

$${\gamma }_2=-{\mathit{pS}}_2\mathcal{F}\text{,}$$

$${\alpha }_1=-2{\mathit{pS}}_4\mathcal{F}\text{,}$$

$${\alpha }_4=\frac{2}{35}\left(7-5S_2-2S_4\right)\mathcal{F}\text{,}$$

$${\alpha }_5=\left[\frac{p}{7}\left(3S_2+4S_4\right)+S_2\right]\mathcal{F}\text{,}$$

$${\alpha }_6=\left[\frac{p}{7}\left(3S_2+4S_4\right)-S_2\right]\mathcal{F}\text{,}$$

where $\mathcal{F}=\frac{k_{\mathrm{B}}T}{D_{\perp }}\rho$, ρ = N/V is the number density of molecules.

On the other hand, a separate calculation of γ₁ requires the averaging of the antisymmetric ST ${\overline{\sigma }}^a=\int d\mathbf{u}f\left(t\text{,}\mathbf{u}\right){\sigma }^a$. This can be done using the Zubarev nonequilibrium statistical operator method [24], and the expression for γ₁ can be written as a certain polynomials of [7], [9], [11], [15]

$${\gamma }_1=\frac{\left(9.54+2.27S_2\right)S_2^2}{2.88+S_2+12.56S_2^2+4.69S_2^3-0.74S_2^4}\mathcal{F}\text{.}$$

The essence of the statistical–mechanical approach which used here for the calculation of γ₁ in uniaxial nematic is that the autocorrelations of the microscopic stress tensor, the stress tensor with the director and fluxes with the order parameter tensor are taken into account.

Having obtained both RVCs γ₁ and γ₂ one can calculate the rest Leslie coefficients as ${\alpha }_2=\frac{1}{2}\left({\gamma }_2-{\gamma }_1\right)$ and ${\alpha }_3=\frac{1}{2}\left({\gamma }_2+{\gamma }_1\right)$.

Thus, according to Eqs. ((2), (3), (4), (5), (6), (7)), γ_i(i = 1, 2) and α_i(i = 1, … , 6) are found to be inversely proportional to the RDC D_⊥, which in turn corresponds to the molecular tumbling in the nematic. We are now in a position to calculate the RVCs γ_i(i = 1, 2) and Leslie coefficients α_i(i = 1, … , 6), provided the temperature dependence of S₂, S₄, and D_⊥, are known. While values of the OPs for cyanobiphenyls are usually fairly easily found [25], the determination of motional constants for rotational diffusion constitutes a formidable task.

Based on the short time expansion of the TCFs [3] ${\mathrm{\Phi }}_{\mathrm{mn}}^L(t)={\mathrm{\Phi }}_{\mathrm{mn}}^L(\infty )+\left[{\mathrm{\Phi }}_{\mathrm{mn}}^L(0)-{\mathrm{\Phi }}_{\mathrm{mn}}^L(\infty )\right]\exp \left(-\frac{t}{{\tau }_{\mathrm{mn}}^L}\right)$, an expression for the orientational correlation time ${\tau }_{\mathrm{mn}}^L$ may be proposed in terms of RDC for a uniaxial molecule. We note that the TCF with m = n = 0 is solely determined by the tumbling motion and ${\tau }_{00}^1$ is related to D_⊥ by

$${\tau }_{00}^1={\left[D_{\perp }\frac{2-2S_2}{1+2S_2}\right]}^{-1}\text{.}$$

Here the first rank (L = 1) functions are relevant for dielectric spectroscopy.

The values of the OPs and the orientational correlation time ${\tau }_{00}^1$ for the hindered rotation of molecules around their molecular short axes can be obtained by means of the MD simulations, based on the realistic atom–atom interaction potentials for 7OCB and HBAB molecules in the nematic phase.

MD simulations for two nematic systems, 7OCB and HBAB, comprising of 128 molecules were carried out at different temperatures, 333, 340, 345 K, for 7OCB, and 345, 353, 360 K, for HBAB, by using the MD program COGNAC [26]. Initially, the LC molecules of each type were placed in a cubic simulation cell with the para axes oriented parallel to the z-axis of the box. The orientations of the molecules were altered in order to prevent a net dipole moment from occurring. During the equilibration procedures the systems were compressed until the experimental densities [27] for 7OCB, and 1 g/cm³ for HBAB, were reached. For each simulation the temperature was kept constant by using a weak coupling to an external heat bath, and the coupling parameter was chosen to be much longer than the integration step and significantly shorter than the simulation length. The NVT equilibrium runs of 5 ns duration were performed using the Berendsen algorithm to control the temperature. After equilibration, the analyses were performed on a trajectory of 2 ns, comprising 5 0000 sets of coordinates. The equations of motion were integrated using the RATTLE algorithm [28] with a time step of 2.0 fs. Simple cutoff method was used in the calculation of electrostatic interaction. Nonbonding interactions were truncated at half of the box length. In the case of 7OCB, the aromatic hydrogens were explicitly considered while the aliphatic CH₂ and CH₃ groups were treated as single interaction centres (united atoms). In the case of HBAB, all CH, CH₂, and CH₃ groups, except for the one between phenyl rings, were treated also as single interaction centres. The Lennard–Jones parameters of two compounds were taken from the optimized potentials for liquid simulations (OPLS). For compounds, 7OCB and HBAB, the force field parameters were taken from Refs.[29], [30], [31]. The fractional atomic charges were determined by using CHELPG [32] scheme following the Hartree–Fock geometry optimization with a 6-31G∗∗ basis set. The charges were calculated with the Gaussian 03 program package [33]. The instantaneous orientation of the director, $\hat{\mathbf{n}}$, was determined from the Cartesian ordering matrix [10] $Q_{\mathrm{zz}}^{\nu \nu }=\frac{1}{N}{\sum }_{\mathrm{j}=\mathrm{i}}\frac{1}{2}\left(3\cos {\theta }_{\mathrm{z}\nu }^j\cos {\theta }_{\mathrm{z}{\nu }^{\prime }}^j-{\delta }_{\nu {\nu }^{\prime }}\right)$, where N is the number of nematic molecules and ${\theta }_{\mathrm{z}\nu }^j$ is the angle between the long molecular axis z and an axis ν fixed in the simulation box. The molecular coordinate system was constructed using eigenvectors of the moment of inertia tensor. The diagonalization of $Q_{\mathrm{zz}}^{\nu \nu }$ gives three eigenvalues; the eigenvector associated with the largest eigenvalue corresponds to the nematic director. For 7OCB and HBAB compounds, the calculated OPs S₂ and S₄, and the experimental values obtained by using the optical technique [25] are given in Table 1(the first three columns). The TCF ${\mathrm{\Phi }}_{00}^1\left(t\right)$ calculated from the trajectory for both compounds, 7OCB and HBAB, at different temperatures (for 7OCB, at T = 333, 340, 345 K, and for HBAB, at T = 345, 353, 360 K) are shown in Figure 1. It should be pointed out that in the infinite limit the TCF is ${\mathrm{\Phi }}_{00}^1\left(\infty \right)=S_2^2$, and one can estimate the average value of the orientational correlation time ${\tau }_{00}^1$. Results for the calculated orientational correlation time ${\tau }_{00}^1$, both for 7OCB and HBAB, and the experimental values for 7OCB, obtained by using the dielectric spectroscopy technique [34] are given in Table 1 (fourth and fifth columns). The RDC D_⊥(T) derived from Eq. (8), for both compounds 7OCB and HBAB, at different temperatures (for 7OCB, at T = 333, 340, 345 K, and for HBAB, at T = 345, 353, 360 K), are shown in Figure 2a. Here T_NI(7OCB) ∼ 348.16 K and T_NI(HBAB) ∼ 373.8 K are the values of both transition temperatures. The temperature dependence of D_⊥(T) can be expressed in an Arrhenius form as$D_{\perp }(T)=D_{\perp }^0\exp \left(-\frac{E}{\mathit{RT}}\right)$, where $D_{\perp }^0$ = 6.09 × 10¹⁰ s⁻¹, and E = 16283.5 J/mol, for 7OCB, and $D_{\perp }^0$ = 0.13 s⁻¹, and E = 48106 J/mol, for HBAB, respectively. Here R is the gas constant.

Having obtained the RDC D_⊥ and the set of OPs S₂ and S₄, one can calculate the rest orientational correlation times, for instance, ${\tau }_{10}^1$ and ${\tau }_{00}^2$ as

$${\tau }_{10}^1={\left[D_{\perp }\frac{2+S_2}{1-S_2}\right]}^{-1}\text{,}$$

$${\tau }_{00}^2={\left[6D_{\perp }\frac{7+5S_2-12S_4}{7+10S_2+18S_4-35S_2^2}\right]}^{-1}\text{.}$$

Results for the orientational correlation times ${\tau }_{10}^1$ and ${\tau }_{00}^2$ are also given in Table 1 (the last two columns).

In order to calculate the values of the RVCs, the length-to-width ratio a = σ_‖/σ_⊥ of both 7OCB and HBAB molecules have been fixed to the value 3.0. The temperature dependencies of the RVCs γ_i(T)(i = 1, 2), calculated using Eqs. (2), (7), for both compounds 7OCB and HBAB, are shown in Figure 2b and c, whereas the temperature dependencies of the Leslie coefficients α_i(i = 1, … , 6), calculated using Eqs. ((3), (4), (5), (6)), are shown in Figure 3, Figure 4, respectively. The experimental values for γ₁ (HBAB) were determined by direct measurements [35], and in good agreement with calculated values, in the temperature range corresponding to the nematic phase of HBAB. Unfortunately, to our knowledge, no systematic measurements of RVC γ₁ or Leslie coefficients α_i(i = 1, … , 6) are reported in the literature for the 7OCB, so that MD predictions can not be reliably tested.

In this letter, we present an investigation of structure, dynamical and viscous properties in the nematic liquid crystals, formed by 4-heptyloxy-4′-cyanobiphenyl (7OCB) and 4-hexyloxy-benzylidene-4′-amino-benzonitrile (HBAB), using a combination of existing statistical–mechanical approaches and molecular dynamics simulations. These simulations has been performed using realistic atom–atom interactions. The orientational order in the simulated liquid crystal compound HBAB ($S_2^{\mathrm{MD}}$ (T = 360 K) = 0.50 ± 0.03) is in good agreement with experimental results derived from the optical method ($S_2^{\mathrm{op}}$ (T = 360 K) = 0.58) [25] and NMR method ($S_2^{\mathrm{NMR}}$ (T = 360 K) = 0.47) [35]. We have determined the values of self-diffusion coefficient corresponding to the hindered rotation both of 7OCB and HBAB molecules around their molecular short axes, as well as the rotational viscosity and Leslie coefficients in the temperature range corresponding to the nematic phase of these compounds. The molecular reorientations in the simulated LCs are in agreement with experimental results derived from broadband spectroscopic methods. In particular, the molecular reorientations were examined using time correlation function, defined by Wigner rotation matrix elements of first rank. From that function we were able to determine the values of the orientational correlation time ${\tau }_{00}^1$ for both compounds 7OCB and HBAB, at a number of temperatures, which can be identified with the experimental dielectric low frequency relaxation time τ_l.f.. Reasonable agreement was found between ${\tau }_{00}^1$ obtained from the TCF and derived from experiments. That time was used for determining of the molecular tumbling, characterized by the rotational self-diffusion coefficient, D_⊥. We have determined the rotational viscosity coefficients γ_i(i = 1, 2) and the Leslie viscosity coefficients α_i(i = 1, … , 6) in the nematic phase. The statistical–mechanical approaches used in calculations of γ_i and α_i rest on a model in which it is assumed that the reorientation of an individual molecule is a stochastic Brownian motion in a certain potential of mean torque. The viscosity coefficients, γ_i and α_i are dependent on the temperature, density and the RSD coefficient, D_⊥. Furthermore, the viscosity is governed by the orientational order in the nematic phase. We found that MD predicts for RVC γ₁ in agreement with experimental values for HBAB. We conclude by pointing out that the combination of theory and computer simulation provides a powerful tool for investigations of viscosity phenomena of real nematics.