## Abstract

The purpose of the work presented in this chapter is to test a recently proven variational principle according to which the irreversible energy dissipation rate is minimal in the linear regime of a nonequilibrium steady state. This test is carried out by performing molecular dynamics simulations of liquid crystals subject to velocity gradients and temperature gradients. Since the energy dissipation rate varies with the orientation of the director of the liquid crystal relative to these gradients and is minimal at certain orientations, this is a stringent test of the variational principle. More particularly, a nematic liquid crystal model based on the Gay-Berne potential, which can be regarded as a Lennard-Jones fluid generalized to elliptical molecular cores, is studied under planar Couette flow, planar elongational flow, and under a temperature gradient. It is found that the director of a nematic liquid crystal consisting of rod-like molecules lies in the vorticity plane at an angle of about 20° to the stream lines in the planar Couette flow. In the elongational flow, it is parallel to the elongation direction, and it is perpendicular to the temperature gradient in a heat flow. These orientations are the ones where the irreversible energy dissipation rate is minimal, so that the variational principle is fulfilled in these three cases.

### Keywords

- liquid crystals
- nonequilibrium molecular dynamics simulation
- shear flow
- elongational flow
- heat conduction
- alignment phenomena
- minimal energy dissipation rate

## 1. Introduction

For a system in thermodynamic equilibrium, there is a variational principle according to which the free energy is minimal, that is, the Helmholtz free energy when the volume, temperature, and number of particles are constant, Gibbs free energy when the pressure, temperature, and number of particles are constant, etc. On the other hand, for systems driven away from equilibrium by an external dissipative field such as a velocity gradient, temperature gradient, chemical potential gradient or an electrical potential gradient doing irreversible work that is converted to heat, there has not been any variational principle to date. However, a theorem originally proposed by Ilya Prigogine stating that a quantity, known as the irreversible energy dissipation rate,

This theorem is not only of basic scientific interest but also of technological and practical interest since shear fields, temperature gradients, concentration gradients, or chemical potential gradients and electrical potential gradients are common examples of external dissipative fields that are ubiquitous in industrial applications and in everyday life. For example, in a lubricated bearing, a planar Couette flow arises in the lubricant in the narrow space between two surfaces rotating at different angular velocities, and

One way of testing this principle is to perform molecular dynamics simulations of microscopic model systems, but then it is hard to find a suitable model system that is easy to analyze. However, liquid crystals are particularly interesting for this purpose because the transport properties and thereby

The simplest kind of liquid crystal is the nematic liquid crystal [3, 4]. It consists of rod-like or plate-like molecules oriented in a certain direction—the director—but there is no translational order, see Figure 1. A nematic liquid crystal cannot support shear stresses, so it is by definition a liquid, but it can support torques, which is the basis for various orientation phenomena relative to external fields. A special case of a nematic liquid crystal is the cholesteric liquid crystal, where the director rotates in space around an axis perpendicular to itself—the cholesteric axis or the optical axis. The spatial rotation period or the pitch is of the order of 1 μm or about 500 molecular diameters. A cholesteric liquid crystal is different from its mirror image, and it is formed by chiral molecules.

There is some theoretical and experimental evidence indicating that the director comes to rest in an orientation where the irreversible energy dissipation rate is minimal in accordance with the variational principle. More specifically, such orientation phenomena have been observed in simulations of shear flow or planar Couette flow [5, 6], in experimental measurements of the viscosity [7] in this flow geometry, and in simulations of planar elongational flow [8]. In the latter case, it is actually possible to prove that the energy dissipation rate must be either minimal or maximal in a steady state in the linear or Newtonian regime by using the linear phenomenological relations between the velocity gradient and the shear stress.

In the case of a nematic liquid crystal subject to a temperature gradient, there are quite a few early experimental works [9, 10, 11, 12, 13, 14] that might imply that the director of a liquid crystal consisting of rod-like molecules orients perpendicularly to this gradient. This means that the heat flow is minimized, since the heat conductivity is minimal in this orientation. Unfortunately, the results of these works are not wholly conclusive because the underlying experiments are very hard to carry out. On the other hand, molecular dynamics simulation of nematic phases of calamitic and discotic soft ellipsoids [15, 16, 17] clearly show that the directors orient perpendicularly and parallel, respectively, to the temperature gradient, so that the heat flow and thereby

This chapter is organized in the following way: in Section 2, the basic theory is outlined, and in Sections 3, 4, and 5, molecular dynamics simulation results and experimental measurements on the director orientation and the irreversible energy dissipation rate are presented and discussed for shear flow or planar Couette flow, planar elongational flow and heat conduction, respectively. In Section 6, the effects of the thermostat are discussed, and finally in Section 7, there is a conclusion. Some background theory is given in the Appendices.

## 2. Basic theory

### 2.1 Order parameter, director, and director angular velocity

In order to describe transport properties of a liquid crystal, we must first define the order parameter, the director, and the director angular velocity. In an axially symmetric system such as a nematic or a smectic *A* liquid crystal, the order parameter, *S*, is given by the largest eigenvalue of the order tensor,

where *N* is the number of particles, and **n**, and it is a measure of the average orientation of the molecules in the system. The order tensor can also be expressed as

The director angular velocity is given by

### 2.2 Director constraint algorithm

Since the molecules studied in the work presented in this chapter are modeled by the Gay-Berne potential, which can be regarded as a Lennard-Jones potential generalized to elliptical molecular cores, see Appendix 2, they are rigid bodies. Therefore, the Euler equations are applied in angular space. Moreover, since the purpose often is to find the stable orientations of the director relative to an external dissipative field, it is interesting to calculate the torque exerted on the liquid crystal, when the director attains various fixed angles relative to this field. This can be done by adding Gaussian constraints to the Euler equations [21],

where *I* is the moment of inertia around the axes perpendicular to the axis of revolution, *i*, *i* by the other molecules, *x*- and *y*-components of the director angular velocity, and *x*- and *y*-components of the director angular acceleration equal to zero. These multipliers are determined in such a way that the director angular acceleration becomes a constant of motion. Then if the initial director angular acceleration and angular velocity are equal to zero, the director will remain fixed in space for all subsequent times and the time averages of the constraint multipliers will be equal to the torque exerted on the director by the external field. Finally, note that the difference between the director angular velocity, **Ω**, and the molecular angular velocities

## 3. Shear flow

### 3.1 The SLLOD equations of motion for shear flow

In order to study shear flow and to calculate the viscosity and director alignment angles relative to the streamlines, it is convenient to apply the SLLOD equations of motion [22]. The name SLLOD stems from the similarity to the Dolls equation of motion derived from the Dolls tensor Hamiltonian. They are synthetic equations of motion that can be used to calculate the viscosity in the linear regime. On the other hand, the idea behind the SLLOD equations of motion is very simple: The velocity of the molecules is divided into the streaming velocity and the thermal velocity. The thermal velocity is related to the temperature, and the streaming velocity is the macroscopic external velocity. The SLLOD equations of motion are an exact description of adiabatic planar Couette flow and a very good approximation of shear flow at constant temperature both in the linear and nonlinear regime. The SLLOD equations are expressed in the following way:

and

where *i*, *m* is the molecular mass, *x*-direction varying linearly in the *z*-direction, see Figure 2, *x*-direction, *i* by the other molecules, and

This expression is obtained by applying Gauss’s principle of least constraint [22]. This principle is essentially the same as the Lagrange’s method for handling constraints. However, Gauss’s principle is more general in that it in addition to constraints involving the molecular coordinates also allows handling of some constraints involving the molecular velocities. This is very useful because it makes it possible to keep the kinetic energy constant whereby the temperature also will be constant. It is possible to show that the ensemble averages of the phase functions and the time correlation functions are essentially canonical when a Gaussian thermostat is applied.

### 3.2 Shear flow of nematic liquid crystals

In a nematic liquid crystal undergoing shear flow, the alignment angle, *θ*, between the director and the streamlines is determined by a mechanical stability criterion, namely, that the antisymmetric pressure must be zero when no external torques act on the system, that is, that the torques exerted by the vorticity and the strain rate cancel out. This makes it possible to derive a relationship between the alignment angle and the viscosity coefficients in the Newtonian regime by using the linear relation between the pressure tensor and the strain rate, see Refs. [3, 4, 23] and Appendix 1,

where

for the preferred alignment angle, *A* phase transition, the ratio

The connection with the variational principle can be made by using the fact that there is an algebraic expression for the irreversible energy dissipation rate,

where the definitions of the viscosity coefficients, *θ* is similar to that given in Figure 3. This function (8) has been obtained by shear flow simulations applying the SLLOD equations of motion [22] to a nematic phase of calamitic soft ellipsoids, see Refs. [5, 6]. The energy dissipation rate is low close to the preferred alignment angle and high when the director is perpendicular to the streamlines and parallel to the velocity gradient. Then, if we study the distribution of the director, we find that it is centered close to the minimum of

or

where

## 4. Planar elongational flow

### 4.1 The SLLOD equations of motion for planar elongational flow

A planar irrotational elongational flow arises when an incompressible liquid expands in the *x*-direction and contracts in the *z*-direction, see Figure 4. Then, the velocity field and the strain rate are given by *x*-direction and contracted in the *z*-direction, the simulation can only continue until the width in the *z*-direction is equal to twice the range of the interaction potential. However, if the angle between the elongation direction and the *x*-axis is set to an angle, *φ*, the periodic lattice of originally quadratic simulation cells is gradually deformed to a lattice of cells shaped like parallelograms. Then, it can be shown that for a special value of this angle, *x*-axis is equal to

and

where *i*, *i* by the other molecules, *m* is the molecular mass, **u** is the streaming velocity, *γ* is the strain rate, **α** is a thermostting multiplier and **β** is a constraint multiplier used to conserve the linear momentum.

### 4.2 Planar elongational flow of nematic liquid crystals

The director alignment angle is in the first place determined by the mechanical stability in the same way as in shear flow whereby the antisymmetric pressure must be zero. In the linear or Newtonian regime, the alignment angle can be found by using the following relation between the antisymmetric pressure and the strain rate, see Appendix 1,

where *θ* denotes the angle between the director and the elongation direction, and

Just as in planar Couette flow or shear flow, the connection to the variational principle can be made by considering the algebraic expression for the irreversible energy dissipation rate in the linear regime,

If the viscosity coefficient *θ* is equal to 45° but this orientation is excluded because of the mechanical stability (12). If *θ* is equal to 0 or 90°, that is, the elongation or contraction direction. Simulations of a nematic phase of calamitic soft ellipsoids have shown that

## 5. Heat conduction

### 5.1 Heat flow algorithm

A temperature gradient can be maintained by keeping different regions, 1 and 2, of a system at different temperatures, see Figure 7. Mathematically, this can be brought about by adding thermostatting terms for each of the regions 1 and 2 to the ordinary Newtonian equations of motion [29],

where *m* is the molecular mass, *i*, and *i* by the other molecules. The thermostatting terms are *α* in Eqs. (4b) and (5), but here, they thermostat the regions 1 and 2 separately. They are determined by applying Gauss’s principle of least constraints using the fact that the weighted kinetic energies are constant:

and

where *ζ* is a multiplier determined in such a way that the linear momentum of the whole system is constant. It goes to zero in the thermodynamic limit.

### 5.2 Heat flow in nematic liquid crystals

The heat flow in an axially symmetric system such as a nematic liquid crystal or a cholesteric liquid crystal is given by

where *T* is the absolute temperature, and **n** is the director. Then, the irreversible energy dissipation rate of the system due to the heat flow becomes,

where the last equality has been obtained by assuming that the director lies in the *zx*-plane forming an angle *θ* with the temperature gradient, see Figure 8. When

The temperature gradient exerts a torque on the molecules around an axis perpendicular to itself and perpendicular to the director, see Figure 8. This torque must be zero in the parallel and perpendicular orientations due to the symmetry, but it is impossible to determine whether these orientations are stable or unstable. Unfortunately, there is no linear relation between the torque and the temperature gradient since they are pseudovectors and polar vectors, respectively, due to the axial symmetry of the system. However, a quantitative relation between them can be obtained by noting that a cross coupling between a pseudo vector and a symmetric second rank tensor is allowed. The latter quantity can be obtained by forming the dyadic product of the temperature gradient, giving the following relation [15],

where **Γ** is the torque density, *μ* is a cross coupling coefficient, and **ε** is the Levi-Civita tensor. The third equality is obtained by assuming that the temperature gradient points in the *z*-direction, and the director lies in the *zx*-plane, see Figure 8, whereby *θ* becomes the angle between these two vectors. This relation fulfills the symmetry condition according to which the torque must be zero when the director is parallel or perpendicular to the temperature gradient. Moreover, the torque is proportional to the square of the temperature gradient for given angle *θ*. Note also that, this relation is analogous to the relation between the strain rate and the antisymmetric pressure in planar elongational flow (12).

The director orientation can be determined by simulating systems, where a temperature gradient and a heat flow are maintained by thermostatting different parts of the system at different temperatures by using the above simulation algorithm (14). Such simulations have shown that the director of nematic liquid crystals consisting of soft calamitic ellipsoids tends to align perpendicularly to the temperature gradient, see Figure 9, whereas the director of nematic liquid crystals consisting of discotic ellipsoids tends to align parallel to the temperature gradient. Thus, the energy dissipation rate is minimal in both cases. Moreover, if the director is constrained to attain a fixed orientation between the perpendicular and parallel orientation relative to temperature gradient by applying the Lagrangian constraint algorithm (3), the torque exerted can be obtained. Then, it is found that this torque turns the director of a calamitic system toward the perpendicular orientation and the director of a discotic system toward the parallel orientation. The same orientation behavior of the directors of calamitic and discotic nematic liquid crystals relative to the temperature gradient was observed in an earlier work [21]. However, then the Evans heat flow algorithm [22] was applied where a fictitious external field under non-Newtonian equations of motion rather than a real temperature gradient drives the heat flow. Therefore, it was not possible to determine whether the orientation phenomena were a real effect or a consequence of the non-Newtonian synthetic equations of motion.

There are also some early experimental works on the orientation of the director of nematic liquid crystals relative to temperature gradients [9, 10, 11, 12, 13, 14] that probably support the conclusions of these heat flow simulations. Unfortunately, it is very difficult to carry out these experiments because if the temperature gradient is too large, there will be convection in the system, and if the temperature gradient is too small, it will not be strong enough to overcome the elastic torques or the surface torques. Therefore, these experiments are not fully conclusive.

Finally, one example where the director orientation relative to a temperature gradient definitely is the one that minimizes the irreversible energy dissipation rate is a cholesteric liquid crystal. In this system, the director rotates in space around the cholesteric axis forming a spiral structure. Then experimental studies, where a temperature gradient is applied, have shown that the cholesteric axis orients parallel to the temperature gradient, whereby the energy dissipation rate is minimized since the heat conductivity is greater in the direction perpendicular to the cholesteric axis than in the parallel direction. Moreover, the whole spiral structure starts rotating in time. This phenomenon is known as thermomechanical coupling [3, 4, 18, 19, 20, 30, 31]. There are quite a few experimental studies available on this phenomenon, where it has been found in a conclusive way that the cholesteric axis remains parallel to the temperature gradient, so this orientation seems to be stable, and thus the irreversible energy dissipation rate is minimal.

We can consequently conclude that the orientation of the director relative to the temperature gradient is consistent with the variational principle [1] even though the coupling between the torque and the temperature gradient is quadratic rather than linear and the system is inhomogeneous. However, the temperature gradient is rather weak, so we still remain in the linear regime.

## 6. Effects of the thermostat

In the above simulations of shear flow and elongational flow, the velocity gradient does work on the system that is converted to heat, which must be removed in order to keep the temperature constant and to maintain a steady state. In a real macroscopic system, this takes place by heat conduction to the system boundaries and this could in principle be arranged in a microscopic simulation cell as well. Unfortunately, this is inconvenient because a temperature gradient of molecular dimensions would make the system inhomogeneous, and thus make it difficult to define the thermodynamic state. Therefore, the temperature is kept constant by forcing the kinetic energy to be a constant of motion by applying a Gaussian thermostat, see Eq. (5). This thermostat was originally devised independently by Hoover *et al*. [32, 33, 34] and by Evans [22]. The equilibrium ensemble averages of the phase functions and time correlation functions generated when this thermostat is applied are essentially canonical [35]. Away from equilibrium, it can be shown that the effect of the Gaussian thermostat on the ensemble averages is proportional to the square of the external field, whereas the thermodynamic fluxes driven by the field are directly proportional to the field in a linear transport process. Thus, the corresponding linear transport coefficients that are equal to the ratio of the flux and the field in the limit of zero field are independent of the thermostat. Therefore, transport coefficients obtained from the simulations of shear flow and elongational flow are independent of the thermostat since there is a linear relation between the velocity gradient and the shear stress in the Newtonian regime and since we are interested in the limit of zero velocity gradient. Neither is the correctness of the variational principle affected by the thermostat since it is valid in the linear regime.

The situation is different in the heat flow simulations because here we actually want a temperature gradient. This gradient is obtained by applying two bar thermostats at different temperatures acting over a limited range and separated by a distance that is long compared to this range, see Figure 7 and Eq. (14). Therefore, the movements of only a small fraction of the molecules are affected by the thermostats, whereas the movements of the majority of the molecules away from the bar thermostats are governed by the ordinary Newtonian equations of motion. Thus, it is reasonable to assume that the influence of the details of the thermostat on the ensemble averages of the phase functions is limited in this case too.

## 7. Conclusion

The purpose of this work has been to test a recently proven variational principle according to which the irreversible energy dissipation rate is minimal in the linear regime of a nonequilibrium steady state. Therefore, we have reviewed molecular dynamics simulations and experimental work on director orientation phenomena in nematic liquid crystals and in cholesteric liquid crystals under external dissipative fields such as velocity gradients and temperature gradients. A general observation that we have made is that in all the examples studied, the director of the liquid crystals seems to attain precisely that alignment angle relative to the external dissipative field that minimizes the irreversible energy dissipation rate.

In a nematic liquid crystal, the director orientation is in the first place determined by a mechanical stability criterion, namely, that the external torques acting on the system must be zero at mechanical equilibrium. This makes it possible to derive an exact relation between the alignment angle relative to the streamlines and the viscosity coefficients in the linear or Newtonian regime of planar elongational flow and of planar Couette flow. Both simulations and experimental measurements imply that the irreversible energy dissipation rate is minimal at this mechanically stable orientation.

It can be shown that the elongation direction is the stable orientation of flow stable calamitic nematic liquid crystals undergoing elongational flow in the linear regime. It can also be shown that the value of the energy dissipation rate is the same in the contraction direction and in the elongation direction, and that this value is either the maximal or the minimal value by using the linear phenomenological relations between the strain rate and the pressure. Simulations of the calamitic soft ellipsoid fluid have shown that the irreversible energy dissipation rate is minimal in the elongation direction.

In calamitic nematic liquid crystals, the heat conductivity is larger in the direction parallel to the director than in the perpendicular direction, and the reverse is true for discotic nematic liquid crystals. Thus, the irreversible energy dissipation rate due to the heat flow depends on the angle between the director and the temperature gradient. When a nematic liquid crystal is subjected to a temperature gradient, a torque is exerted on the molecules. Due to symmetry, this torque must be proportional to the square of the temperature gradient and it must be zero when the director is parallel or perpendicular to this gradient.

In simulations of nematic phases of soft ellipsoids under a temperature gradient, it turns out that the director of a calamitic nematic liquid crystal aligns perpendicularly to the temperature gradient, whereas the director of a discotic nematic liquid crystal attains the parallel orientation. In both cases, the irreversible energy dissipation rate is minimal. These simulation results are probably supported by some experimental measurements, but they are difficult to carry out in practice so they are not fully conclusive.

Finally, one system where there is definitely a conclusive experimental evidence for the fact that the director attains the orientation that minimizes the energy dissipation rate due to a temperature gradient is the cholesteric liquid crystal. The cholesteric axis of droplets of cholesteric liquid crystals orient parallel to a temperature gradient and the director rotates. This is a well-established phenomenon observed in studies of thermomechanical coupling, and since the heat conductivity is lower in the direction of the cholesteric axis than in the perpendicular direction, the energy dissipation rate is minimal in this case.

Thus, the director orientation relative to a temperature gradient also follows the variational principle even though there is a quadratic coupling between the torque and the temperature gradient. However, the temperature gradients are rather low so we are still in the linear regime.

## Acknowledgments

We gratefully acknowledge financial support from the Knut and Alice Wallenberg Foundation (Project number KAW 2012.0078) and Vetenskapsrådet (Swedish Research Council) (Project number 2013-5171). The simulations were performed using resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC, HPC2N, and NSC.) We also acknowledge PRACE for awarding us access to Hazelhen based in Germany at Rechenzentrum Stuttgart.

The relation between the velocity gradient, **P**, is more complicated in an axially symmetric system such as nematic liquid crystal than in an isotropic fluid due to the lower symmetry. In order to derive the linear phenomenological relations between the velocity gradient and the pressure, it is appropriate to begin by identifying the thermodynamic forces and fluxes in the expression for the irreversible entropy production [3, 4, 23, 36]:

where *T* is the absolute temperature, and **u** is the streaming velocity. The various parts of the second rank tensor are denoted in the following manner: the symmetric traceless part is given by **ε** is the Levi-Civita tensor. Three pairs of thermodynamic forces and fluxes can be identified by inspection of the irreversible entropy production, namely, the symmetric traceless pressure tensor and the traceless strain rate,

and

where the products involving the Levi-Civita tensor **ε** are defined in the following way: *ζ*, which, according to the ORR, is equal to the cross coupling coefficient *κ* between the traceless strain rate and the difference between the trace of the pressure tensor and the equilibrium pressure.

Application of a planar Couette velocity gradient, *zx*-plane at an angle *θ* relative to the stream lines, see Figure 2, by application of an electric or a magnetic field gives the following relations between the pressure tensor components and the strain rate in a director based coordinate system

and

where

In a planar elongational flow [8, 26, 27, 28], where the elongation direction is parallel to the *x*-axis, the contraction direction is parallel to the negative *z*-axis, and the velocity field is equal to *θ* is the angle between the director and the elongation direction or *x*-axis,

and

If these expressions for the pressure tensor are inserted in the expression for energy dissipation rate (A.1), we obtain

for planar Couette flow and

for planar elongational flow. The subscript *γ* denotes that the average is evaluated in a nonequilibrium ensemble at a finite shear rate.

In order to evaluate the above expressions for the irreversible work in shear flow, elongational flow, and heat flow, we have simulated a coarse grained model system composed of molecules interacting via a purely repulsive version of the commonly used Gay-Berne potential [16, 17, 21],

where

and

where the parameter

The denominators in Eqs. (A.8a) and (A.8b) are never equal to zero because the absolute value of the scalar product