Variational Principle for Nonequilibrium Steady States Tested by Molecular Dynamics Simulation of Model Liquid Crystal Systems

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.

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, ${\dot{w}}_{\mathit{irr}}$ , is minimal in the linear regime of a nonequilibrium steady state has recently been proven [1]. This quantity is defined as the irreversible work per unit time and unit volume that is done by a dissipative external field on the system [2]. Thus, there is a variational principle for nonequilibrium steady states.

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 ${\dot{w}}_{\mathit{irr}}$ is equal to the product of the shear rate and the shear stress. Another example is the heat flow between a hot region and a cold region such as the inside and the outside of a building. Then, ${\dot{w}}_{\mathit{irr}}$ is equal to the product of the heat flow and the temperature gradient. Still another example is an electric heating element where an electric potential difference or voltage drives an electric current, and ${\dot{w}}_{\mathit{irr}}$ is equal to the product of the voltage and the current. Finally, chemical potential gradients arise when various substances are mixed and they begin to diffuse, and ${\dot{w}}_{\mathit{irr}}$ is equal to the product of the chemical potential gradients and the matter currents.

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 ${\dot{w}}_{\mathit{irr}}$ depend on the orientation relative to the streamlines or the temperature gradient, and at certain orientations, ${\dot{w}}_{\mathit{irr}}$ is minimal. Thus, it can be determined whether these orientations actually are attained by the liquid crystal. Moreover, it is possible to orient the liquid crystal in an experimental measurement by applying an electric or magnetic field and in molecular dynamics simulations by applying a constraint torque. This means that ${\dot{w}}_{\mathit{irr}}$ can be measured or calculated as a function of the orientation relative to the dissipative field.

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.

Figure 1.

A nematic phase of the Gay-Berne fluid undergoing planar Couette flow. The velocity gradient is directed in the vertical direction and the streamlines are directed in the horizontal direction. Note that the director forms an angle with streamlines of about 18°.

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 ${\dot{w}}_{\mathit{irr}}$ are minimized. However, one system, where the director definitely orients perpendicularly to the temperature gradient, is the cholesteric liquid crystal, where the cholesteric axis orients parallel to the temperature gradient, so that the director becomes perpendicular to this gradient, and the heat flow is minimized [3, 4, 18, 19, 20], which is in agreement with the variational principle.

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,

$$\mathbf{Q}=\frac{3}{2}\left(\frac{1}{N}\sum _{i=1}^N{\widehat{\mathbf{u}}}_i{\widehat{\mathbf{u}}}_i-\frac{1}{3}\mathbf{1}\right),$$ E1

where N is the number of particles, and $\left\{{\widehat{\mathbf{u}}}_i;1\le i\le N\right\}$ is some characteristic vector of the molecule; in the case of bodies of revolution, it can be taken to be parallel to the axis of revolution but in a more realistic all atom model some other vector in the molecule has to be defined as ${\widehat{\mathbf{u}}}_i$ , and $\mathbf{1}$ is the unit second rank tensor. When the molecules are perfectly aligned in the same direction, the order parameter is equal to unity, and when the orientation is completely random, it is equal to zero. The eigenvector corresponding to the order parameter is defined as the director, n, and it is a measure of the average orientation of the molecules in the system. The order tensor can also be expressed as

$$\mathbf{Q}=\frac{3}{2}S\left(\mathbf{nn}-\frac{1}{3}\mathbf{1}\right).$$ E2

The director angular velocity is given by $\mathbf{\Omega }=\mathbf{n}\times \dot{\mathbf{n}}$ . In a macroscopic system, the order tensor and the order parameter are functions of the position in space, but in a small system such as a simulation cell with dimensions of the order of some ten molecular lengths, there is only one director and one order parameter for the whole system.

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],

$$I{\dot{\mathbf{\omega }}}_i={\mathbf{\Gamma }}_i+{\lambda }_x\frac{\partial {\mathrm{\Omega }}_x}{\partial {\mathbf{\omega }}_i}+{\lambda }_y\frac{\partial {\mathrm{\Omega }}_y}{\partial {\mathbf{\omega }}_i},$$ E3

where I is the moment of inertia around the axes perpendicular to the axis of revolution, ${\mathbf{\omega }}_i$ is the angular velocity of molecule i, ${\mathbf{\Gamma }}_i$ is the torque exerted on molecule i by the other molecules, ${\mathrm{\Omega }}_x$ and ${\mathrm{\Omega }}_y$ are the x- and y-components of the director angular velocity, and ${\lambda }_x$ and ${\lambda }_y$ are Gaussian constraint multipliers keeping the 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 $\left\{{\mathbf{\omega }}_i;1\le i\le N\right\}$ ; the director angular velocity can be regarded as the angular velocity of the average orientation of the molecules. If the director angular velocity is constrained to be zero by applying Eq. (3), the molecular angular velocities are still nonzero and the right hand side of Eq. (3) is nonzero. The Gaussian constraint simply forces the molecules to rotate in such a way that average orientation stays the same.

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:

$${\dot{\mathbf{r}}}_i=\frac{{\mathbf{p}}_i}{m}+\gamma r_{\mathit{zi}}{\mathbf{e}}_x$$ E4a

and

$${\dot{\mathbf{p}}}_i={\mathbf{F}}_i-\gamma p_{\mathit{zi}}{\mathbf{e}}_x-\alpha {\mathbf{p}}_i,$$ E4b

where ${\mathbf{r}}_i$ and ${\mathbf{p}}_i$ are the position and peculiar momentum, that is, the momentum relative to the streaming velocity, of molecule i, m is the molecular mass, $\gamma =\partial u_x/\partial z$ is the shear rate, that is, there is a streaming velocity $u_x$ in the x-direction varying linearly in the z-direction, see Figure 2, ${\mathbf{e}}_x$ is the unit vector in the x-direction, ${\mathbf{F}}_i$ is the force exerted on molecule i by the other molecules, and $\alpha$ is a thermostatting multiplier given by the constraint that the linear peculiar kinetic energy is a constant of motion,

$$\alpha =\frac{{\sum }_{i=1}^N\left[{\mathbf{F}}_i\cdot {\mathbf{p}}_i-\gamma p_{\mathit{ix}}{p}_{\mathit{iz}}\right]}{{\sum }_{i=1}^N{\mathbf{p}}_i^2}.$$ E5

Figure 2.

Planar Couette flow or shear flow arises when there is a streaming velocity in the x-direction, varying linearly in the z-direction, $\mathbf{u}=\mathit{\gamma z}{\mathbf{e}}_x$ , where $\gamma =\partial u_x/\partial z$ is the shear rate or velocity gradient. The expression for the relation between the velocity gradient and the pressure tensor becomes simpler by using a director-based coordinate system $\left({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3\right)$ , where the director n points in the ${\mathbf{e}}_3$ -direction, obtained by rotating the ordinary laboratory-based coordinate system $\left({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3\right)$ with an angle θ around the ${\mathbf{e}}_y={\mathbf{e}}_2$ -axis. Reproduced from Ref. [6] with the permission of AIP Publishing.

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,

$$⟨p_2^a⟩=-{\tilde{\gamma }}_1\frac{\gamma }{4}-{\tilde{\gamma }}_2\frac{\gamma }{4}cos2\theta =0,$$ E6

where ${\tilde{\gamma }}_1$ is the twist viscosity, ${\tilde{\gamma }}_2$ is the cross coupling coefficient between the antisymmetric pressure and the strain rate, and $⟨p_2^a⟩$ is the antisymmetric pressure in the vorticity direction perpendicular to the streamlines and perpendicular to the velocity gradient. The angular brackets denote that the pressure tensor is the ensemble average of a phase function. Then, if $⟨p_2^a⟩$ is equal to zero, we obtain

$$cos2{\theta }_0=-{\tilde{\gamma }}_1/{\tilde{\gamma }}_2,$$ E7

for the preferred alignment angle, ${\theta }_0$ , provided that the ratio $\left|{\tilde{\gamma }}_1/{\tilde{\gamma }}_2\right|$ is less than one. Then the liquid crystal is said to be flow stable. This condition is fulfilled in many liquid crystals, and ${\theta }_0$ is between 10 and 20° both in real systems and in simplified coarse grained model systems such as the soft ellipsoid fluid, see Refs. [5, 6, 7], Figure 1, and Appendix 2. Note, however, that for some systems, often near the nematic-smectic A phase transition, the ratio $\left|{\tilde{\gamma }}_1/{\tilde{\gamma }}_2\right|$ is greater than one. This means that there is no orientation angle where the antisymmetric pressure is zero, so that no steady state is attained. Then the liquid crystal is said to be flow unstable and the director will rotate forever [3, 4, 23, 24].

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, ${\dot{w}}_{\mathit{irr}}$ , of a flow stable nematic liquid crystal, given by the dyadic product of the symmetric traceless pressure, $\overline{\mathbf{P}}$ , and the traceless strain rate, $\overline{\mathbf{\nabla u}}$ ,

$${\dot{w}}_{\mathit{irr}}=-\overline{\mathbf{P}}:\overline{\mathbf{\nabla u}}=\left(\eta +\frac{{\tilde{\eta }}_1}{6}+\frac{{\tilde{\eta }}_3}{2}{sin}^{2}2\theta +\frac{{\tilde{\eta }}_2}{2}cos2\theta \right){\gamma }^2,$$ E8

where the definitions of the viscosity coefficients, $\eta$ , ${\tilde{\eta }}_1$ , ${\tilde{\eta }}_2$ , and ${\tilde{\eta }}_3$ _, and the derivation are given in Appendix 1. If the values of the various viscosity coefficients are inserted, it is found that the functional dependence of ${\dot{w}}_{\mathit{irr}}$ on θ 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 ${\dot{w}}_{\mathit{irr}}$ . This minimum has also been observed in simulations of shearing nematic phases of discotic soft ellipsoids [5] and when experimentally measured, viscosity coefficients are inserted in the Eqs. (7) and (8) and the resulting alignment angle is determined [7]. Thus the system seems to select the alignment angle that minimizes the irreversible energy dissipation rate in accordance with the variational principle. This also means that the energy dissipation rate (8) must be minimal at the preferred alignment angle, ${\theta }_0$ , given by Eq. (7). Thus, the derivative of the function (8) with respect to $\theta$ must be zero for ${\theta }_0$ , giving an additional relation between the viscosity coefficients and the alignment angle,

$$cos2{\theta }_0-=\frac{{\tilde{\eta }}_2}{2{\tilde{\eta }}_3}$$ E9

or

$$2{\tilde{\eta }}_3{\tilde{\gamma }}_1+{\tilde{\eta }}_2{\tilde{\gamma }}_2=0,$$ E10

where ${\theta }_0$ has been eliminated by using Eq. (7). The Eqs. (7) and (9) do not coincide but they must still give the same value of ${\theta }_0$ . This provides an important cross check for the correctness of the simulation algorithms and experimental methods used to determine the viscosity coefficients and for the computer programs used to run the simulations.

Figure 3.

The irreversible energy dissipation rate, ${\dot{w}}_{\mathit{irr}}$ , Eq. (8), due to the strain rate of a nematic liquid crystal phase of calamitic soft ellipsoids as a function of the director alignment angle, θ, is obtained by using the Gaussian constraint algorithm (3) to fix the director at various angles relative to the streamlines. The preferred alignment angle attained when no constraints are applied is equal to about 20° which is close to the minimum of ${\dot{w}}_{\mathit{irr}}$ . Reproduced from Ref. [6] with the permission of AIP Publishing.

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 $\mathbf{u}=\gamma \left(x{\mathbf{e}}_x-z{\mathbf{e}}_z\right)$ and $\nabla \mathbf{u}=\overline{\nabla \mathbf{u}}=\gamma \left({\mathbf{e}}_x{\mathbf{e}}_x-{\mathbf{e}}_z{\mathbf{e}}_z\right)$ . Planar elongational flow can be studied by applying a special version of the SLLOD equations of motion. Then the problem is that, if the simulation cell is elongated in the 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, ${\phi }_0$ , the lattice of parallelograms can be remapped onto the original quadratic lattice after a certain time period, so that the simulation becomes continuous, that is, the Kraynik-Reinelt boundary conditions, see Figure 5 and Refs. [8, 25, 26, 27, 28] for details. Then, if the angle between the elongation direction and the x-axis is equal to ${\phi }_0$ , the velocity gradient becomes $\mathbf{\nabla u}=\gamma \left({\mathbf{e}}_x^{\prime }{\mathbf{e}}_x^{\prime }-{\mathbf{e}}_z^{\prime }{\mathbf{e}}_z^{\prime }\right)$ , where ${\mathbf{e}}_x^{\prime }={\mathbf{e}}_{x}cos{\phi }_0-{\mathbf{e}}_{z}sin{\phi }_0$ and ${\mathbf{e}}_z^{\prime }={\mathbf{e}}_{x}sin{\phi }_0+{\mathbf{e}}_{y}cos{\phi }_0$ are the elongation and contraction directions. Inserting this gradient in the SLLOD equations of motion gives,

$${\dot{\mathbf{r}}}_i=\frac{{\mathbf{p}}_i}{m}+{\mathbf{r}}_i\cdot \mathbf{\nabla u}=\frac{{\mathbf{p}}_i}{m}+\gamma {\mathbf{r}}_i\cdot \left({\mathbf{e}}_x^{\prime }{\mathbf{e}}_x^{\prime }-{\mathbf{e}}_z^{\prime }{\mathbf{e}}_z^{\prime }\right)$$ E11a

and

$${\dot{\mathbf{p}}}_i={\mathbf{F}}_i-{\mathbf{p}}_i\cdot \mathbf{\nabla u}-\alpha {\mathbf{p}}_i-\mathbf{\beta }={\mathbf{F}}_i-\gamma {\mathbf{p}}_i\cdot \left({\mathbf{e}}_x^{\prime }{\mathbf{e}}_x^{\prime }-{\mathbf{e}}_z^{\prime }{\mathbf{e}}_z^{\prime }\right)-\alpha {\mathbf{p}}_i-\mathbf{\beta },$$ E11b

where ${\mathbf{r}}_i$ and ${\mathbf{p}}_i$ are the position and peculiar momentum, that is, the momentum relative to the macroscopic streaming velocity, of molecule i, ${\mathbf{F}}_i$ is the force exerted on molecule 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.

Figure 4.

Schematic representation of a nematic phase of a soft ellipsoid fluid undergoing irrotational extensional flow. The system is elongated in the horizontal direction and contracted in the vertical direction while the volume is constant. The molecules tend to be aligned in the elongation direction.

Figure 5.

The Kraynik-Reinelt boundary conditions. The original simulation cell is square 1. When the angle between the elongation direction and the horizontal direction is equal to ${\phi }_0$ , square 1 is deformed to a parallelogram, which, after a given time interval, becomes the dashed parallelogram, partially covering the squares (1–6). Then the triangles a’, b’, and c’ in the parallelogram are periodic copies of the triangles a, b, and c in square 1. If the primed triangles are moved to the corresponding unprimed triangles, a square is recovered and the simulation can continue. Reproduced from Ref. [8], with permission from the PCCP Owner Societies.

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,

$$⟨p_2^a⟩=-{\tilde{\gamma }}_2\frac{\gamma }{2}sin2\theta ,$$ E12

where θ denotes the angle between the director and the elongation direction, and ${\tilde{\gamma }}_2$ is the cross coupling coefficient between the antisymmetric pressure and the strain rate. From this expression, it is apparent that the alignment angle must be either 0 or 90°, that is, where the torque exerted by the strain rate is equal to zero. For a flow stable calamitic nematic liquid crystal, the cross coupling coefficient ${\tilde{\gamma }}_2$ is negative [6], so that the 0° orientation parallel to the elongation direction is mechanically stable, and the 90° orientation is unstable.

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,

$${\dot{w}}_{\mathit{irr}}=-\overline{\mathbf{P}}:\overline{\mathbf{\nabla u}}=\left(4\eta +\frac{2{\tilde{\eta }}_1}{3}+2{\tilde{\eta }}_3{cos}^{2}2\theta \right){\gamma }^2.$$ E13

If the viscosity coefficient ${\tilde{\eta }}_3$ is positive, this expression is minimal when θ is equal to 45° but this orientation is excluded because of the mechanical stability (12). If ${\tilde{\eta }}_3$ on the other hand is negative, this expression attains the same minimal value when θ 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 ${\tilde{\eta }}_3$ is less than zero [8], so that the energy dissipation rate is minimal in the stable orientation also in this case of planar elongational flow. This is in agreement with the variational principle [1]. See also Figure 6 where the angular distribution of the director around the elongation direction is displayed.

Figure 6.

The angular distribution, $p\left(\theta \right)$ , of the director of a calamitic nematic liquid crystal consisting of soft ellipsoids around the elongation direction where the angle between the director and the elongation direction is denoted by θ. Reproduced from Ref. [8], with permission from the PCCP Owner Societies.

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],

$$m{\ddot{\mathbf{r}}}_i={\mathbf{F}}_i-{\widehat{w}}_{1i}{\alpha }_{1}m{\dot{\mathbf{r}}}_i-{\widehat{w}}_{2i}{\alpha }_{2}m{\dot{\mathbf{r}}}_i-\zeta ,$$ E14

where m is the molecular mass, ${\dot{\mathbf{r}}}_i$ and ${\ddot{\mathbf{r}}}_i$ are the velocity and the acceleration of molecule i, and ${\mathbf{F}}_i$ is the force exerted on molecule i by the other molecules. The thermostatting terms are ${\widehat{w}}_{1i}{\alpha }_{1}m{\dot{\mathbf{r}}}_i$ and ${\widehat{w}}_{2i}{\alpha }_{2}m{\dot{\mathbf{r}}}_i$ where ${\widehat{w}}_{1i}$ and ${\widehat{w}}_{2i}$ are two normalized weight functions. These terms are actually similar to the thermostatting term in Eq. (4b), but here region 1 and region 2 are thermostatted separately at different temperatures. This is achieved by letting the weight functions be Gaussian functions centered in region 1 and 2, respectively, and with decay lengths that are considerably shorter than the distance between these two regions. In this way, only the molecules in region 1 contribute to the temperature in region 1, and only the molecules in region 2 contribute to the temperature in region 2. The molecules far away from the centers of these regions move according to the ordinary Newtonian equations of motion. Note that, it is not necessary to use Gaussian weight functions; it is possible to use any function with a maximum and a rather short decay length. The parameters ${\alpha }_1$ and ${\alpha }_2$ are thermostatting multipliers in the same way as the multiplier α 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:

$$\frac{1}{2}\sum _{i=1}^N{\widehat{w}}_{1i}m{\dot{\mathbf{r}}}_i^2=E_{k1}$$ E15a

and

$$\frac{1}{2}\sum _{i=1}^N{\widehat{w}}_{2i}m{\dot{\mathbf{r}}}_i^2=E_{k2},$$ E15b

where $E_{k1}$ and $E_{k2}$ are the weighted kinetic energies for region 1 and 2, respectively. The algebraic expressions for the thermostatting multipliers are given in Ref. [15]. The parameter ζ 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.

Figure 7.

A temperature gradient is maintained by thermostatting one region (dark gray) of the system at a high temperature and another region (light gray) at a low temperature, whereby heat will flow from the high temperature region to the low temperature region. Reproduced from Ref. [15], with permission from the PCCP Owner Societies.

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

$$⟨{\mathbf{J}}_Q⟩=-\left[{\lambda }_{\Vert \Vert }\mathbf{nn}+{\lambda }_{\perp \perp }\left(\mathbf{1}-\mathbf{nn}\right)\right]\cdot \frac{\mathbf{\nabla }T}{T},$$ E16

where $⟨{\mathbf{J}}_Q⟩$ is the heat current density, ${\lambda }_{\Vert \Vert }$ is the heat conductivity parallel to the director of an ordinary achiral nematic liquid crystal or parallel to the cholesteric axis of a cholesteric liquid crystal, ${\lambda }_{\perp \perp }$ is the heat conductivity perpendicular to the director of a nematic liquid crystal or perpendicular to the cholesteric axis of a cholesteric liquid crystal, 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,

$$\begin{array}{c}{\dot{w}}_{\mathit{irr}}=-⟨{\mathbf{J}}_Q⟩\cdot \frac{\mathbf{\nabla }T}{T}=\frac{1}{T^2}\left[{\lambda }_{\perp \perp }\mathbf{\nabla }T\cdot \mathbf{\nabla }T+\left({\lambda }_{\Vert \Vert }-{\lambda }_{\perp \perp }\right){\left(\mathbf{n}\cdot \mathbf{\nabla }T\right)}^2\right]\\ ={\left|\frac{{\partial }_{z}T}{T}\right|}^2\left[{\lambda }_{\perp \perp }+\left({\lambda }_{\Vert \Vert }-{\lambda }_{\perp \perp }\right){cos}^2\theta \right],\end{array}$$ E17

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 ${\lambda }_{\Vert \Vert }>{\lambda }_{\perp \perp }$ , as in a nematic liquid crystal consisting of calamitic molecules, the heat current density and thereby ${\dot{w}}_{\mathit{irr}}$ are maximal when the temperature gradient and the director are parallel and minimal when they are perpendicular to each other. Conversely, when ${\lambda }_{\Vert \Vert }<{\lambda }_{\perp \perp }$ , as in a nematic liquid crystal consisting of discotic molecules, the heat current density and the irreversible energy dissipation rate are maximal when the director is perpendicular to the temperature gradient and minimal when it is parallel to the temperature gradient.

Figure 8.

A schematic view of a nematic liquid crystal subject to a temperature gradient is shown. The temperature gradient ∇T points in the z-direction, and the director n lies in the zx-plane forming an angle θ to the z-axis. Then a torque Γ arises in the direction of the y-axis. Reproduced from Ref. [15], with permission from the PCCP Owner Societies.

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],

$$\mathbf{\Gamma }=\mu \mathbf{\epsilon }:\mathbf{nn}\cdot \frac{\nabla T}{T}\frac{\nabla T}{T}=\mu \mathbf{n}\cdot \frac{\nabla T}{T}\mathbf{n}\times \frac{\nabla T}{T}=\mu {\left|\frac{{\partial }_{z}T}{T}\right|}^{2}cos\theta sin\theta {\mathbf{e}}_y=\frac{1}{2}\mu {\left|\frac{{\partial }_{z}T}{T}\right|}^{2}sin2\theta {\mathbf{e}}_y,$$ E18

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.

Figure 9.

The angular distribution, $p\left(\theta \right)$ , of the director of a nematic liquid crystal consisting of soft calamitic ellipsoids around the temperature gradient where the angle between the director and the temperature gradient is denoted by θ. Reproduced from Ref. [15], with permission from the PCCP Owner Societies.

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, $\mathbf{\nabla u}$ , and the pressure tensor, 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]:

$$\sigma =-\frac{1}{T}\left\{\overline{\mathbf{P}}:\overline{\mathbf{\nabla u}}+2{\mathbf{P}}^a\mathbf{\cdot }\left(\mathrm{½}\mathbf{\nabla }\times \mathbf{u}-\mathbf{\Omega }\right)+\left(\frac{1}{3}\mathrm{Tr}\left(\mathbf{P}\right)-p_{\mathit{eq}}\right)\mathbf{\nabla }\mathbf{\cdot }\mathbf{u}\right\},$$ EA.1

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 $\overline{\mathbf{A}}=\mathrm{½}\left(\mathbf{A}+{\mathbf{A}}^T\right)-\left(1/3\right)\mathrm{Tr}\left(\mathbf{A}\right)\mathbf{1}$ and the pseudovector dual of the antisymmetric part is denoted by ${\mathbf{A}}^a=-\mathrm{½}\mathbf{\epsilon }\mathbf{:}\mathbf{A}=-\mathrm{½}{\epsilon }_{\mathit{\alpha \beta \gamma }}{\mathrm{A}}_{\mathit{\gamma \beta }}$ , where ε 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, $\overline{\mathbf{P}}$ and $\overline{\mathbf{\nabla u}}$ , the antisymmetric pressure and the difference between the rotation and the director angular velocity, ${\mathbf{P}}^a$ and $\mathrm{½}\mathbf{\nabla }\times \mathbf{u}-\mathbf{\Omega }$ , and the difference between the trace of the pressure tensor and the equilibrium pressure of a quiescent liquid crystal, and the trace of the strain rate, $\left(1/3\right)\mathrm{Tr}\left(\mathbf{P}\right)-p_{\mathit{eq}}$ and $\mathbf{\nabla }\mathbf{\cdot }\mathbf{u}$ . Note that the strain rate is defined as $\mathrm{½}\left[\mathbf{\nabla u}+{\left(\mathbf{\nabla u}\right)}^T\right]$ , and it is always symmetric. In a uniaxially symmetric nematic liquid crystal, the relations between the pressure and the velocity gradient can be deduced by symmetry arguments, and they can be expressed in a few different equivalent ways [23, 36]. It has been found that a notation due to Hess [36] is the most convenient one for deducing Green-Kubo relations and NEMD-algorithms:

$$⟨\overline{\mathbf{P}}⟩=-2\eta \overline{\nabla \mathbf{u}}-{\tilde{\eta }}_1\overline{\overline{\mathbf{nn}}\cdot \overline{\nabla \mathbf{u}}}-2{\tilde{\eta }}_3\overline{\mathbf{nn}}\overline{\mathbf{nn}}\mathbf{:}\overline{\nabla \mathbf{u}}+2{\tilde{\eta }}_2\overline{\overline{\mathbf{nn}}\cdot \mathbf{\epsilon }\cdot \left(\mathrm{½}\nabla \times \mathbf{u}-\mathbf{\Omega }\right)}-\zeta \overline{\mathbf{nn}}\nabla \cdot \mathbf{u},$$ EA.2a

$$⟨{\mathbf{P}}^a⟩=-\frac{{\tilde{\gamma }}_1}{2}\left(\mathrm{½}\mathbf{\nabla }\times \mathbf{u}-\mathbf{\Omega }\right)-\frac{{\tilde{\gamma }}_2}{2}\mathbf{\epsilon }:\left(\overline{\mathbf{nn}}\cdot \overline{\nabla \mathbf{u}}\right),$$ EA.2b

and

$$\frac{1}{3}⟨\mathrm{Tr}\left(\mathbf{P}\right)⟩-p_{\mathit{eq}}=-{\eta }_V\mathbf{\nabla }\cdot \mathbf{u}-\kappa \overline{\mathbf{nn}}\mathbf{:}\overline{\nabla \mathbf{u}},$$ EA.2c

where the products involving the Levi-Civita tensor ε are defined in the following way: $\mathbf{\epsilon }:\mathbf{A}={\epsilon }_{\mathit{\alpha \beta \gamma }}{\mathrm{A}}_{\mathit{\gamma \beta }}$ and $\mathbf{A}\cdot \mathbf{\epsilon }\cdot \mathbf{B}={\mathrm{A}}_{\mathit{\alpha \beta }}{\epsilon }_{\mathit{\beta \gamma \delta }}{B}_{\delta }$ . The quantities $\eta$ , ${\tilde{\eta }}_1$ , and ${\tilde{\eta }}_3$ are shear viscosities, ${\tilde{\gamma }}_1$ is the twist viscosity, ${\eta }_V$ is the volume viscosity, ${\tilde{\eta }}_2$ is the cross coupling coefficient relating the difference between the rotation and the director angular velocity, and the symmetric traceless pressure. According to the Onsager reciprocity relations (ORR), this coefficient is equal to ${\tilde{\gamma }}_2/2$ , the cross coupling coefficient relating the traceless strain rate and the antisymmetric pressure. The trace of the strain rate and the symmetric traceless pressure are related by the cross coupling coefficient, ζ, 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, $\mathbf{\nabla u}=\gamma {\mathbf{e}}_z{\mathbf{e}}_x$ , where $\gamma ={\partial }_z{u}_x$ is the shear rate and fixation of the director in the 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 $\left({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3\right)$ where the director points in the ${\mathbf{e}}_3$ -direction:

$$⟨{\overline{p}}_{11}⟩=\left(\eta +\frac{{\tilde{\eta }}_3}{3}\right)\gamma sin2\theta ,$$ EA.3a

$$⟨{\overline{p}}_{22}⟩=\frac{1}{3}\left({\tilde{\eta }}_1+{\tilde{\eta }}_3\right)\gamma sin2\theta ,$$ EA.3b

$$⟨{\overline{p}}_{33}⟩=-\left(\eta +\frac{{\tilde{\eta }}_1}{3}+\frac{2{\tilde{\eta }}_3}{3}\right)\gamma sin2\theta ,$$ EA.3c

$$⟨{\overline{p}}_{31}⟩=\left(\eta +\frac{{\tilde{\eta }}_1}{6}\right)\gamma cos2\theta +{\tilde{\eta }}_2\frac{\gamma }{2},$$ EA.3d

and

$$2⟨p_2^a⟩=⟨{\widehat{\lambda }}_2⟩=-{\tilde{\gamma }}_1\frac{\gamma }{2}-{\tilde{\gamma }}_2\frac{\gamma }{2}cos2\theta ,$$ EA.3e

where ${\widehat{\lambda }}_2$ is the external torque density acting on the system. From these equations, it is apparent that the various elements of the pressure tensor are linear functions of $sin2\theta$ and $cos2\theta$ , so the various viscosity coefficients can be evaluated by fixing the director at a few different angles relative to the stream lines and calculating the averages of the pressure tensor elements.

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 $\mathbf{u}=\gamma \left(x{\mathbf{e}}_x-z{\mathbf{e}}_z\right)$ , so that the velocity gradient becomes $\mathbf{\nabla u}=\gamma \left({\mathbf{e}}_x{\mathbf{e}}_x-{\mathbf{e}}_z{\mathbf{e}}_z\right)$ . Then the linear relations between the velocity gradient and the pressure become the following in a director-based coordinate system $\left({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3\right)$ where the director points in the ${\mathbf{e}}_1$ -direction, and θ is the angle between the director and the elongation direction or x-axis, ${\mathbf{e}}_2={\mathbf{e}}_y$ and ${\mathbf{e}}_3={\mathbf{e}}_1\times {\mathbf{e}}_2$ ,

$$⟨{\overline{p}}_{11}⟩=-2\left(\eta +\frac{{\tilde{\eta }}_3}{3}\right)\gamma cos2\theta ,$$ EA.4a

$$⟨{\overline{p}}_{22}⟩=-\frac{2}{3}\left({\tilde{\eta }}_1+{\tilde{\eta }}_3\right)\gamma cos2\theta ,$$ EA.4b

$$⟨{\overline{p}}_{33}⟩=2\left(\eta +\frac{{\tilde{\eta }}_1}{3}+\frac{2{\tilde{\eta }}_3}{3}\right)\gamma cos2\theta ,$$ EA.4c

$$⟨{\overline{p}}_{31}⟩=\left(2\eta +\frac{{\tilde{\eta }}_1}{3}\right)\gamma sin2\theta ,$$ EA.4d

and

$$2⟨p_2^a⟩=⟨{\widehat{\lambda }}_2⟩=-{\tilde{\gamma }}_2\gamma sin2\theta .$$ EA.4e

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

$$\begin{array}{c}{\dot{w}}_{\mathit{irr}}=\overline{\mathbf{P}}:\overline{\mathbf{\nabla u}}=-{⟨\frac{1}{2}\left({\overline{p}}_{33}^s-{\overline{p}}_{11}^s\right)sin2\theta -{\overline{p}}_{31}^{s}cos2\theta ⟩}_{\gamma }\\ =\left(\eta +\frac{{\tilde{\eta }}_1}{6}+\frac{{\tilde{\eta }}_3}{2}{sin}^{2}2\theta +\frac{{\tilde{\eta }}_2}{2}cos2\theta \right){\gamma }^2,\end{array}$$ EA.5

for planar Couette flow and

$${\dot{w}}_{\mathit{irr}}=\left(4\eta +\frac{2{\tilde{\eta }}_1}{3}+2{\tilde{\eta }}_3{sin}^{2}2\theta \right){\gamma }^2$$ EA.6

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],

$$U\left({\mathbf{r}}_{12},{\widehat{\mathbf{u}}}_1,{\widehat{\mathbf{u}}}_2\right)=4\epsilon \left({\widehat{\mathbf{r}}}_{12},{\widehat{\mathbf{u}}}_1,{\widehat{\mathbf{u}}}_2\right){\left(\frac{{\sigma }_0}{r_{12}-\sigma \left({\widehat{\mathbf{r}}}_{12},{\widehat{\mathbf{u}}}_1,{\widehat{\mathbf{u}}}_2\right)+{\sigma }_0}\right)}^{18},$$ EA.7

where ${\mathbf{r}}_{12}={\mathbf{r}}_2-{\mathbf{r}}_1$ is the distance vector from the center of mass of molecule 1 to the center of mass of molecule 2, ${\widehat{\mathbf{r}}}_{12}$ is the unit vector in the direction of ${\mathbf{r}}_{12}$ , $r_{12}$ is the length of the vector ${\mathbf{r}}_{12}$ _, and ${\widehat{\mathbf{u}}}_1$ and ${\widehat{\mathbf{u}}}_2$ are the unit vectors parallel to the axes of revolution of molecule 1 and molecule 2. The parameter ${\sigma }_0$ is the length of the axis perpendicular to the axis of revolution, that is, the minor axis of a calamitic ellipsoid of revolution and the major axis of a discotic ellipsoid of revolution. The strength and range parameters are given by

$$\epsilon \left({\widehat{\mathbf{r}}}_{12},{\widehat{\mathbf{u}}}_1,{\widehat{\mathbf{u}}}_2\right)={\epsilon }_0{\left[1-{\chi }^2{\left({\widehat{\mathbf{u}}}_1\cdot {\widehat{\mathbf{u}}}_2\right)}^2\right]}^{-1/2}{\left\{1-\frac{{\chi }^{\prime }}{2}\left[\frac{{\left({\widehat{\mathbf{r}}}_{12}\cdot {\widehat{\mathbf{u}}}_1+{\widehat{\mathbf{r}}}_{12}\cdot {\widehat{\mathbf{u}}}_2\right)}^2}{1+{\chi }^{\prime }{\widehat{\mathbf{u}}}_1\cdot {\widehat{\mathbf{u}}}_2}+\frac{{\left({\widehat{\mathbf{r}}}_{12}\cdot {\widehat{\mathbf{u}}}_1-{\widehat{\mathbf{r}}}_{12}\cdot {\widehat{\mathbf{u}}}_2\right)}^2}{1-{\chi }^{\prime }{\widehat{\mathbf{u}}}_1\cdot {\widehat{\mathbf{u}}}_2}\right]\right\}}^2$$ EA.8a

and

$$\sigma \left({\widehat{\mathbf{r}}}_{12},{\widehat{\mathbf{u}}}_1,{\widehat{\mathbf{u}}}_2\right)={\sigma }_0{\left\{1-\frac{\chi }{2}\left[\frac{{\left({\widehat{\mathbf{r}}}_{12}\cdot {\widehat{\mathbf{u}}}_1+{\widehat{\mathbf{r}}}_{12}\cdot {\widehat{\mathbf{u}}}_2\right)}^2}{1+\chi {\widehat{\mathbf{u}}}_1\cdot {\widehat{\mathbf{u}}}_2}+\frac{{\left({\widehat{\mathbf{r}}}_{12}\cdot {\widehat{\mathbf{u}}}_1-{\widehat{\mathbf{r}}}_{12}\cdot {\widehat{\mathbf{u}}}_2\right)}^2}{1-\chi {\widehat{\mathbf{u}}}_1\cdot {\widehat{\mathbf{u}}}_2}\right]\right\}}^{-1/2},$$ EA.8b

where the parameter $\chi$ is equal to $\left({\kappa }^2-1\right)/\left({\kappa }^2+1\right)$ and $\kappa$ is the ratio of the axis of revolution and the axis perpendicular to this axis, ${\chi }^{\prime }$ is equal to $\left({{\kappa }^{\prime }}^{1/2}-1\right)/\left({{\kappa }^{\prime }}^{1/2}+1\right)$ and ${\kappa }^{\prime }$ is the ratio of the potential energy minima of the side by side and end to end configurations of calamitic ellipsoids or the ratio of the edge-to-edge and face-to-face configurations of discotic ellipsoids, and ${\epsilon }_0$ denotes the depth of the potential minimum in the cross configuration, where ${\widehat{\mathbf{r}}}_{12}$ , ${\widehat{\mathbf{u}}}_1$ _, and ${\widehat{\mathbf{u}}}_2$ are perpendicular to each other. The parameters $\kappa$ and ${\kappa }^{\prime }$ have been given the values 3 and 5, respectively, for the calamitic ellipsoids and 1/3 and 1/5 for the discotic ellipsoids.

The denominators in Eqs. (A.8a) and (A.8b) are never equal to zero because the absolute value of the scalar product ${\widehat{\mathbf{u}}}_1\cdot {\widehat{\mathbf{u}}}_2$ is less than or equal to one since ${\widehat{\mathbf{u}}}_1$ and ${\widehat{\mathbf{u}}}_2$ are unit vectors, and the absolute values of the parameters $\chi$ and ${\chi }^{\prime }$ are less than one. The ordinary Lennard-Jones potential is recovered in the limit when $\kappa$ and ${\kappa }^{\prime }$ go to one. Note that, the potential is purely repulsive, so there are no potential minima but the value of ${\kappa }^{\prime }$ optimized for the attractive Gay-Berne potential has been retained. The transport properties of this system of purely repulsive soft ellipsoids are similar to those of a system where the molecules interact according to the conventional Gay-Berne potential with attraction as well, so the results are still relevant.