By means of a combination of equilibrium Monte Carlo and molecular dynamics simulations and nonequilibrium molecular dynamics we investigate the ordered, uniaxial phases (i.e., nematic and smectic A) of a model liquid crystal. We characterize equilibrium behavior through their diffusive behavior and elastic properties. As one approaches the equilibrium isotropic-nematic phase transition, diffusion becomes anisotropic in that self-diffusion *D*_{⊥} in the direction orthogonal to a molecule’s long axis is more hindered than self-diffusion *D*_{∥} in the direction parallel to that axis. Close to nematic-smectic A phase transition the opposite is true, *D*_{∥} < *D*_{⊥}. The Frank elastic constants *K*_{1}, *K*_{2}, and *K*_{3} for the respective splay, twist, and bend deformations of the director field $n\u0302$ are no longer equal and exhibit a temperature dependence observed experimentally for cyanobiphenyls. Under nonequilibrium conditions, a pressure gradient applied to the smectic A phase generates Poiseuille-like or plug flow depending on whether the convective velocity is parallel or orthogonal to the plane of smectic layers. We find that in Poiseuille-like flow the viscosity of the smectic A phase is higher than in plug flow. This can be rationalized via the velocity-field component in the direction of the flow. In a sufficiently strong flow these smectic layers are not destroyed but significantly bent.

## I. INTRODUCTION

Liquid crystals (LCs) have been a source of fascinatingly abundant physics for more than a century. It all stems from a simple molecular fact: the rigidity of their molecules gives rise to anisotropic interactions, that is, their potential energy depends on their mutual orientations. Macroscopically, this anisotropy leads to a series of phase transitions that break rotational and translational symmetries in a step-wise fashion. For example, as the temperature is lowered from an isotropic LC below the nematic transition, a macroscopic molecular orientation emerges. At an even lower temperature, one-dimensional layering of the molecular centers of mass emerges, namely the smectic A phase.

Because of their capacity to reorient, also in response to external fields, LCs are used in a wide range of applications: from electronic displays, to microlasers^{1–3} and lubricants.^{4} LCs are ubiquitous in biological settings. Nucleic acids, proteins, carbohydrates, and fats exhibit liquid-crystalline mesophases.^{5,6} Biological LCs possess both structural and functional properties. Biological membranes, for example, are smectic A analogues^{7,8} and have flexoelectric responses,^{6,9} thus coupling sensing and actuation.^{6} Cytoskeletal filaments can experience spatial order and alignment both at the level of the mesh size (≃10 nm) and of the whole cell (≃10 *μ*m), leading to short- and long-range directionality.^{10,11} Finally, LCs have also risen to an important role in biomedical sciences and applications,^{12,13} and in our comprehension of morphogenesis and evolution of living organisms.^{14}

Hydrodynamic flow of an LC fundamentally perturbs its equilibrium properties which is manifested in the interplay of preferential orientation, surface anchoring, topological defects, and flow velocity. Considerable interest in the flow of LCs has grown^{15} for at least two reasons: (i) it gives access to fundamental quantities as the viscosity coefficients and (ii) novel methods of manipulation of liquids at the microscale and lab-on-a-chip technologies can be ushered by the anisotropy of LC interactions.

Although the hydrodynamics of nematic LCs has been well studied,^{16,17} surprisingly much less is known about the next LC phase with lower symmetry, that is, the smectic A phase. Among the first experimental investigations of the smectic A rheological properties is the work of Porter *et al.* in 1966.^{18} Helfrich introduced the concept of permeation.^{19} When a smectic A phase flows in a capillary with velocity perpendicular to the layers, the molecules will ‘permeate’ through the immobile smectic pattern, producing a plug flow. The fundamental hydrodynamic equations for smectics were written down by Martin, Parodi, and Pershan^{20} and later considered in simplified form by de Gennes.^{21} More recently, Bennett and Hess^{22} studied the Miesowicz, Helfrich, and Leslie viscosities in proximity of the nematic to smectic phase transition. Walker and Stewart^{23} introduced a new theoretical approach to smectic A flow. In addition, shear flow experiments on systems close to the nematic-smectic transition or well into the smectic phases have found complex transformations of the layer structure. Typically, different orientations of the layers coexist or the layer normal is reoriented from the direction of the shear gradient to the vorticity direction.^{24–27} Alternatively, undulations of the layers are observed.^{25}

The experimental picture up to date is broad and well understood on a macroscopic length scale.^{27–29} Additionally, theoretical approaches such as elasticity theory capture the experimental picture adequately.^{21,30} Moreover, the interest in equilibrium^{31,32} and nonequilibrium^{33,34} properties of LC confined to mesoscopic or even nanoscopic geometries has grown in the past years.^{35} On this length scale the elasticity theory does not provide a correct quantitative description anymore, such that molecular simulations are a vital tool to further explore this emerging field.

Thus, there is a great need to have model potentials for LCs that are both sufficiently realistic and computationally not too demanding. For example, the well-established Gay-Berne model of an LC (studied, for example, in Ref. 22) is based essentially on a Lennard-Jones potential where the orientation dependence of the interaction is built-in via an orientation dependent van der Waals radius and depth of the attractive well. It is this latter feature which renders the Gay-Berne model particularly cumbersome. This is a consequence of the relatively large shape anisotropy of the molecules which are a notorious problem in computer simulations. In fact, not only the computational effort of the rather complex potential but also long equilibration phases limit the size of simulations. Because we are interested in simulation setups that mimic experimental conditions, we are faced with the challenge to devise a model potential for a LC that reproduces its phenomenology but that at the same time allows to study large systems. We therefore introduce a novel model for smectic A fluids (Section II). In Section III we describe the details of the Monte Carlo and molecular dynamics simulations that we have performed. Section IV is given to a description of equilibrium and nonequilibrium properties of our model smectic fluid. Finally, we discuss our results and conclusions in Section V.

## II. MODEL

### A. Mesogens

We consider a LC composed of *N* molecules (i.e., mesogens) interacting via a pairwise additive potential such that the total configurational potential energy may be cast as

where *r*_{ij} = *r*_{i} − *r*_{j} is the distance vector connecting the centers of mass of mesogens *i* and *j* located at *r*_{i} and *r*_{j}, respectively. In addition to *r*_{ij}, the interaction potential *φ*_{mm} depends on the orientations *ω*_{ i} and *ω _{j}* of the interacting mesogenic pair. In Eq. (2.1) we use shorthand notation $R=r1,\u2026,rN$ and $\Omega =\omega 1,\u2026,\omega N$ to indicate that Φ

_{mm}depends on the configuration of the

*N*mesogens. On account of the uniaxial symmetry of the mesogens, $\omega =\theta ,\varphi $ where

*θ*and

*ϕ*are the polar and azimuthal angle specifying the orientation of the mesogens’ symmetry axes in a space-fixed frame of reference.

Because *φ*_{mm} is orientation dependent, we split it into an isotropic and an anisotropic contribution according to

where $rij=rij$. More specifically, we take

where ε is the depth of the attractive well, *σ* is the van der Waals diameter of the spherically symmetric mesogenic core, *λ* is the inverse Debye screening length of the attractive Yukawa tail, and subscripts “rep” and “att” refer to repulsive and attractive contributions to the isotropic-core potential, respectively. We emphasize that Eq. (2.3) differs from *φ*_{iso} used in previous work.^{36,37} In fact, here we take a Yukawa-like instead of the Lennard-Jones potential for reasons to be discussed later. Parameters

are introduced to guarantee that the minimum of *φ*_{iso} remains at $rijmin=\sigma $ and that $\phi iso(rijmin)=\epsilon /10$ regardless of *λ*.

For the anisotropic interactions, we follow Giura and Schoen^{38} and expand *φ*_{anis} in terms of the so-called rotational invariants.^{39} To proceed we employ symmetry considerations, namely the invariance of *φ*_{anis} with respect to the transformations $\omega \u2009i\u2192\omega \u2009i\u2032=\u2212\omega \u2009i$ or $\omega j\u2192\omega j\u2032=\u2212\omega j$ as well as its invariance upon replacing *r*_{ij} by −*r*_{ij}. One may then cast *φ*_{anis} as

where $r\u0302ij=rij/rij$. Throughout this work the caret is used to indicate a unit vector. If one limits the expansion of *φ*_{anis} to the three leading terms,^{38} the anisotropy function can be expressed as

where $P2x=123x2\u22121$ is the second Legendre polynomial and the (dimensionless) anisotropy parameters 2ε_{1} = − ε_{2} = 4.0 are fixed throughout this work. Last but not least, we take *λσ* = 3.0 such that *φ*_{iso} is short-range.

At this stage, a couple of comments are worthwhile. The term proportional to ε_{1} in Eq. (2.6) corresponds to the orientation dependence of interactions in the celebrated Maier–Saupe model.^{16} Notice that for fixed orientations of mesogens *i* and *j* and for ε_{2} = 0, *φ*_{mm} is isotropic. Nevertheless, nematic phases may form because configurations in which $u\u0302(\omega \u2009i)\u22c5u\u0302(\omega j)=\xb11$ (head-tail symmetry) are energetically most favorable.^{38} However, just including the Maier–Saupe term in Eq. (2.6) would be insufficient if the formation of smectic A phases is desired. In order to support both long-range orientational and one-dimensional positional order, the term proportional to ε_{2} is indispensable. More specifically, terms proportional to ε_{2} can be understood as ellipsoidal eccentricity of the mesogens’ shape. Notice that the anisotropy parameters taken here are larger than in previous work.^{36,37} In order to balance these rather strong anisotropic interactions, we take advantage of the Yukawa-like potential [see Eq. (2.3)] and use a weak [i.e., $\phi iso(rijmin)=\epsilon /10$] and short-ranged (i.e., *λσ* = 3.0) potential. One immediately realizes that for our present choice of ε_{1} > 0 and ε_{2} < 0 together with $\epsilon 1<\epsilon 2$, *φ*_{anis} is attractive on the one hand as long as $u\u0302(\omega \u2009i)\u22c5u\u0302(\omega j)=\xb11$ *and* in addition $u\u0302(\omega \u2009i)\u22c5r\u0302ij=u\u0302(\omega j)\u22c5r\u0302ij=0$. If, on the other hand, $u\u0302(\omega \u2009i)\u22c5u\u0302(\omega j)=u\u0302(\omega \u2009i)\u22c5r\u0302ij=u\u0302(\omega j)\u22c5r\u0302ij=\xb11$, *φ*_{anis} is repulsive and therefore layer formation as in a smectic A structure is supported. These features are illustrated by the plot in Figure 1.

### B. Confining walls

To investigate the impact of flow on ordered liquid-crystalline structures, it is necessary to control the direction in which this order emerges with respect to the direction of flow. Experimentally, this is commonly achieved by placing the liquid crystal between solid substrates. The surfaces of these substrates are prepared in such a way that they anchor the mesogens in their vicinity in a certain way. “Anchoring” refers to an energetic discrimination of desirable orientations of the mesogens with respect to the plane of the substrate. The solid surfaces can either be prepared by the deposition of chemicals^{40,41} or mechanically by rubbing the solid surface in particular ways.^{42–45}

We adopt this experimental concept here and place our liquid crystal between plane parallel, atomically corrugated solid substrates that are separated by a distance *s*_{z} along the *z*-axis. The contribution of the mesogen-substrate interaction to the total configurational potential energy can then be cast as

assuming each substrate to be composed of a monolayer of *N*_{s} rigidly fixed atoms. In each monolayer, the atoms are arranged according to the (100) structure of the face-centered cubic lattice. We take the lattice constant to be given by $\u2113/\sigma =43$ corresponding to an areal density of the substrate atoms of *ρ*_{s}*σ*^{2} = 2/ℓ^{2}. Taking the configuration of substrate atoms ** S** = (

*s*_{1}, …,

*s*_{2Ns}) to be fixed, $rij\u2032=ri\u2212sj$ is the distance between the center of mass of mesogen

*i*and substrate atom

*j*. Moreover, the substrates are in registry, that is, each atom in one substrate is exactly opposite to its counterpart in the other substrate. In Eq. (2.7),

where the distance dependence of the interactions is described by the well-known Lennard–Jones potential.

The expression in Eq. (2.8) differs from the standard Lennard–Jones potential by the anchoring function 0 ≤ *g*_{α} ≤ 1 involved in the latter. One may view the anchoring function as a mathematical “device” to realize alignment of the mesogens with specific directions. In the smectic A phase, this allows us to control the formation of individual layers such that the layer normal points in a desired direction with respect to the direction of flow. Throughout this work we will always keep the direction of flow along the *x*-axis and therefore parallel to both substrates. Here, we consider

where $e\u0302\alpha $ (*α* = x or y) is the versor (i.e., a unit vector) in the *α*-direction. From Eqs. (2.9) and (2.8), it is immediately clear that mesogen-substrate attraction is fully “switched on” if $u\u0302$ is parallel to the versor; on the contrary, if $u\u0302$ is orthogonal to the versor, fluid-substrate attraction is completely suppressed.

Finally, before concluding this section a few comments seem appropriate. First, the atomically corrugated wall employed here is computationally more demanding than a smooth one. This is because in the latter case *φ*_{ms} turns out to be a one-body potential whereas the computation of Φ_{ms} from Eq. (2.7) involves at most the evaluation of a double sum over 2*NN*_{s} terms. However, atomic corrugation is indispensable because under flow, mesogens sufficiently close to the substrate surfaces need to experience friction for the entire system to eventually reach a steady state.

Second, an inspection of Figure 1 and of Eq. (2.8) reveals that for *g*_{α} = 1 the attractive wells of *ϕ*_{mm} and *ϕ*_{ms} are equally deep as far as side-side configurations of the mesogens are concerned. This is crucial because then the mesogen-substrate potential is on the one hand sufficiently strong to control the global nematic director $n\u03020$ over the course of a simulation, but, on the other hand, sufficiently weak to prevent the system (or a part of it) from becoming dynamically arrested, as we have verified.

## III. COMPUTATIONAL DETAILS

We begin our presentation of results by first discussing bulk properties of the (three-dimensional) liquid-crystal model introduced in Section II A. These results are obtained from MC simulations as far as static properties are concerned and from equilibrium molecular dynamics (EMD) simulations where we concentrate on self-diffusion as the key dynamic property at equilibrium. Both types of simulations are carried out in the isothermal-isobaric ensemble.

For the MC simulations we employ an adapted Metropolis algorithm as, for instance, described in the book by Allen and Tildesley.^{46} On account of the anisotropic shape of the mesogens, the generation of a Markov chain proceeds in two consecutive steps. First, it is decided with equal probability whether to randomly displace the center of mass of a sequentially chosen mesogen by a small amount or to rotate it by a small angle increment about a randomly chosen axis. Both processes are realized through a Metropolis algorithm.^{46}

Once all *N* mesogens have been considered, one attempt is made to change the volume *V* of the system by a small amount. Again, an energy criterion is employed to decide whether to accept or reject that volume change.^{46} The ratio *N* : 1 between attempted displacements/rotations and volume change is prompted by the fact that the latter is computationally much more demanding. This is because for *V* → *V*′, $R\u2192R\u2032=RV/V\u20323$ which then requires a recalculation of all *N* (*N* − 1)/2 terms in Eq. (2.1) *in principle*. For displacements and rotations, this number is reduced to just *N* as one can easily verify.

In practice, however, both these numbers are reduced in *absolute* but not in *relative* terms because we employ a potential cutoff *r*_{c} = 3.0*σ* beyond which we neglect *φ*_{mm}; no correction is applied for those neglected terms. Whereas this already reduces the computational burden substantially, we also employ a combination of a link-cell and Verlet neighborlist.^{46} Thus, a mesogen is considered a neighbor of another one if their center-of-mass distance is smaller than or equal to *r*_{N} = 3.8*σ*. This reduces the computational cost associated with the search for mesogens that are within the radius *r*_{c} of the cutoff sphere centered on a reference mesogen.

The sequence of *N* + 1 attempts to change the position/orientation of mesogens and system volume constitutes a MC cycle. Our results are based upon 5 × 10^{3} such cycles for equilibration followed by another 5 × 10^{4} cycles during which ensemble averages are being computed. These runs are sufficiently long because they have always been started from quite well equilibrated configurations as we have verified in a few cases.

In the complementary EMD simulations, we integrate the classical equations of motion numerically employing the velocity Verlet algorithm proposed by Ilnytskyi and Wilson.^{47} Forces and torques required by this algorithm are computed by differentiating Φ_{mm} in Eq. (2.1) with respect to *r*_{i} and *ω _{i}*, respectively. To generate trajectories in phase space that are compatible with the isothermal-isobaric ensemble, we employ a combination of a Nosé-Hoover thermostat

^{48–50}and a Hoover barostat.

^{46}Similar to the MC simulations, a potential cutoff and a neighborlist are employed here where the same values of

*r*

_{c}and

*r*

_{N}are used. We equilibrate the system for 2.0 × 10

^{5}to 3.0 × 10

^{5}time steps followed by another 10

^{6}such steps during which time averages are collected. In the equilibrium simulations (MC and EMD) our systems comprise between 500 and 5000 mesogens.

After characterizing various mesophases of our LC model through its equilibrium properties in Section IV A, we then turn to a discussion of its nonequilibrium properties. We focus here on the smectic A phase that we establish below for our model at sufficiently low *T*. Our NEMD simulations are carried out exclusively with *N* = 2.0 × 10^{4} mesogens placed in a rectangular nanofluidic channel of volume *V* = *s*_{x}*s*_{y}*s*_{z}, where all three side lengths *s*_{α} ≈ 24. Along the *z*-axis the liquid crystal is confined by solid substrates where the mesogen-substrate interaction potential is given in Eq. (2.8).

Based upon equilibrium results discussed in Section IV A we choose *T* = 0.70 and *P* = 1.00 in the isothermal-isobaric NEMD simulations to make sure that the liquid crystal is sufficiently deep in the smectic A phase (see also Figure 2). Under these conditions we obtain a mean number density of *ρσ*^{3} ≃ 1.40.

In the NEMD simulations the equations of motion are integrated numerically using again the velocity Verlet algorithm proposed by Ilnytskyi and Wilson^{47} with a time step of *δt* = 10^{−3}. A steady-state flow is initiated and maintained by applying a body force $Fe=Fee\u0302x$ to each mesogen and at every time step. Every simulation is carried out according to the following protocol:

We equilibrate the system in the isothermal-isobaric ensemble in the absence of flow.

A hydrodynamic flow is applied in a second equilibration run under isochoric rather than isobaric conditions until a steady state has been reached.

As in the second equilibration run, observables of interest are finally recorded in a production run under isochoric conditions.

Because of the external body force, the system needs to be thermostatted permanently to achieve a stationary nonequilibrium state. To that end, we employ a local version of the Nosé-Hoover thermostat.^{48–50} To implement this thermostat, one first needs to realize that the symmetry of the system is broken along the *z*-axis on account of the solid substrates. Therefore, we discretize the system along the *z*-axis into slices of thickness 2.4*σ* and apply the thermostat to each slice separately and independently. This assures that momentum is conserved locally and no viscous heating occurs throughout the system.

Finally, because we are not concerned with any specific material all quantities will be given in dimensionsless (i.e., “reduced”) units. For example, length will be expressed in units of *σ*, energy in units of ε, mass in units of *m*, and temperature in units of ε/*k*_{B}, where *k*_{B} is Boltzmann’s constant. Other quantities such as time, pressure, or density can be expressed in terms of combinations of the basic ones.^{46}

## IV. RESULTS

### A. Phases at thermodynamic equilibrium

As the focus of this study is on ordered (i.e., nematic and smectic A) phases, we introduce measures of order first. To that end we follow Eppenga and Frenkel^{51} and introduce the alignment tensor

in its local version, where *δ* denotes the Dirac *δ*-function and

is the local density. In Eq. (4.1), **Q** is a traceless, symmetric second-rank tensor where ⊗ denotes the tensor (i.e., dyadic) product, **1** is the unit tensor, and $\u2026$ indicates an ensemble average in the isothermal isobaric ensemble.

By analogy with Eq. (4.1) we also introduce the nonlocal (volume averaged) tensor

which satisfies the eigenvalue equation

and where $n\u0302k$ is the *k*th eigenvector corresponding to eigenvalue *λ _{k}*.

Because **Q**_{V} can be represented by a 3 × 3 matrix, the secular equation corresponding to Eq. (4.4) is a cubic polynomial in *λ _{k}*. Its roots can be found analytically using Cardano’s

^{52}formula.

^{53}In an infinitely large system the discriminant

*D*

_{3}of the cubic secular equation vanishes indicating that all of its roots are real and two of them are equal. In any system of finite size,

*D*

_{3}< 0 and therefore the three eigenvalues are different. One can interpret this as an ostensible biaxality which strongly depends on system size.

^{37}

Solving instead Eq. (4.4) with the ensemble averaged tensor $QV\u2032=QV$ gives the ensemble averaged eigenvalues $\lambda k\u2032$, where in general $\lambda k=\lambda k\u2032$. Notice, however, that on account of the head-tail symmetry of the mesogens, $n\u0302k\u2032\u2260n\u0302k$. Following Eppenga and Frenkel we then define the nematic order parameter via $S=maxk=13\lambda k\u2032$ and take the associated eigenvector as the nematic director $n\u03020$.

In the smectic A phase mesogens align themselves on average with $n\u03020$. In addition, one observes the formation of layers where the layer normal also points along $n\u03020$. Therefore, to quantify the degree of layer formation we introduce the smectic order parameter via

This form of Λ is inspired by perceiving the local density in a smectic A phase as a plane wave evolving in the direction of $n\u03020$. In Eq. (4.5), the *a priori* unknown parameter *d* is a measure of layer spacing. It is closely related to the length of a mesogen’s long axis (see Figure 1). However, the definition of the latter is somewhat ambiguous in our model. This is because *φ*_{mm} has no zeros for end-end configurations of a pair of mesogens. Moreover, on account of thermal fluctuations the layer spacing may vary slightly across the smectic A phase. Hence, we adjust *d* for every configuration $R,\Omega $ such that Λ is maximized.^{54} Over the range of temperatures for which the smectic A phase is thermodynamically stable we obtain *d* ≃ 1.58.

The temperature dependence of *S* and Λ across the isotropic-nematic (IN) and the nematic-smectic A (NSmA) phase transitions is illustrated by the plots shown in Figure 2. Starting at high *T*, both *S* and Λ nearly vanish. They are not exactly zero on account of a finite-size effect in the isotropic phase. In the case of *S* this finite-size effect has been rationalized and analyzed in depth by Greschek and Schoen^{37} who investigated a model system closely related to the present one.

A similar finite size effect is anticipated for the smectic order parameter Λ. Consider first a smectic A phase at *T* = 0 assuming perfect smectic layers and no thermal fluctuations. In this case, $rj\u22c5n\u03020/d=n$, where *n* ∈ ℕ. Hence, the real part of the complex exponential in Eq. (4.5) is one whereas the imaginary part vanishes. On the contrary, $rj\u22c5n\u03020/d$ is not an integer but a real number in the absence of smectic layers (for example, in the nematic phase where $n\u03020$ is still well defined). If in this case the limit *N* → ∞ is considered, $\u2211j=1\u221e\u2192\u222b01dxexp2\pi ix=0$. Hence for finite *N*, for which the sum in Eq. (4.5) cannot be approximated by such an integral, a small residual value of Λ remains.

As *T* is lowered, *S* suddenly rises indicating the IN phase transition. Notice that for those *T* for which *S* is already quite substantial, Λ ≤ 0.10. Hence, even though orientational order is already long-range, positional order is not. The smectic order parameter begins to rise at even lower *T* for which *S* already begins to level off.

To locate the temperature *T*_{IN} of the IN phase transition more precisely, we follow Greschek and Schoen’s work and perform a finite-size scaling analysis.^{37} To that end we solve Eq. (4.4). Using the middle eigenvalue *λ*_{0} we then define the second-order Binder cumulant as

Greschek and Schoen^{37} found and subsequently rationalized that for a sufficiently weak discontinuous phase transition $g20\u22121$ for different *N* intersect at a unique value of the thermodynamic field (*P* or *T*) driving the transition.

Plots in Figure 3 show that in the isotropic phase cumulants for different system sizes are spread out. In fact, here $g20\u221dN$.^{37} Thus, the larger *N* the larger is $g20\u22121$. In the nematic phase this dependence on system size is inverted. As one can see from Figure 3, a unique intersection of all curves exists demarcating the temperature *T*_{IN} ≃ 0.87 of the IN phase transition.

The plots in Figure 3 also indicate that order-parameter fluctuations are decreasing substantially as one penetrates deeper into the nematic phase. This is because $g20\u22121$ is a measure of these fluctuations *relative* to the mean value of the order parameter as one can easily verify from Eq. (4.6). Over the temperature range considered in Figure 3, these fluctuations decay by about four orders of magnitude between the isotropic and the nematic phase.

Moreover, *T*_{IN} from the cumulant analysis agrees quite well with the estimated inflection point of the curve $ST$ plotted in Figure 2. At *T*_{IN}, *S*_{IN} ≃ 0.35 from MC which agrees very nicely with the universal value of $SIN=13$ predicted by Landau-de Gennes theory.^{55} Adopting the inflection point as an alternative, operational definition of the value of *T* at which a phase transition takes place, we obtain from the data for Λ plotted in Figure 2, *T*_{NSmA} ≃ 0.78 for the temperature at the NSmA phase transition.

To further characterize disordered and ordered liquid-crystalline phases in our model, we now turn to a discussion of self-diffusion. Here, it seems interesting to distinguish between self-diffusion in the direction to a mesogen’s long axis and self-diffusion perpendicular to that axis. Let us therefore define $ri\u2225=ri\u22c5u\u0302\omega iu\u0302\omega i$ and $ri\u22a5=ri\u2212ri\u2225$. With these quantities we introduce the mean-square displacements (MSDs),

where $\u2026t0$ indicates an average over a set of different initial times $t0$ due to stationarity^{56} of the self-diffusion process. With Eq. (4.7) we define

as quantitative measures of self-diffusion at thermodynamic equilibrium.

Plots of $\Delta r\u22a52$ in Figure 4(a) indicate that regardless of the nature of the phase in question, the MSD exhibits two different regimes. At short times all three MSDs exhibit ballistic motion of the mesogens, that is, $\Delta r\u22a52\u221dt2$. Ballistic motion arises as long as intermolecular interactions do not matter. Consequently, all three data sets plotted in Figure 4(a) can be represented by a master curve for *t* ≲ 0.1.

The time range over which the mesogens move ballistically depends on the density of the system and the kinetic energy of the mesogens. Hence, somewhere in the range 0.1 ≲ *t* ≲ 1.0 intermolecular interactions become important and the motion of mesogens changes from ballistic to diffusive. The latter is characterized by $\Delta r\u22a52\u221dt$ which then persists as *t* → ∞.

From plots in Figure 4(b), one notices similar features of the MSD at short times and in the long-time regime as far as isotropic and nematic phases are concerned. However, for the smectic A phase the MSD changes from ballistic to subdiffusive first (i.e., $\Delta r\u22252\u221dt\nu $, *ν* < 1) before eventually becoming diffusive again at long times. The subdiffusive regime at intermediate times reflects the one-dimensional positional ordering in the smectic A phase. In other words, self-diffusion in the direction of $n\u03020$ is hindered because of the one-dimensional solid-like structure in that direction. However, because the structure within each smectic layer remains fluidic, diffusive behavior eventually sets in but only at markedly longer times compared with both the isotropic and the nematic phase.

Because all MSDs plotted in Figure 4 exhibit diffusive behavior at sufficiently long times, we compute *D*_{∥} and *D*_{⊥} from Eq. (4.8). Plots in Figure 5(a) show that with decreasing *T*, *D*_{∥} and *D*_{⊥} are monotonically decreasing functions. In the isotropic phase, *D*_{∥} ≈ *D*_{⊥} which makes sense because no spatial direction is distinguished here.

As one enters the nematic phase, one sees from Figure 5(b) that at higher *T* diffusional isotropy is no longer preserved. More specifically, diffusion along the mesogens’ axes is larger than the diffusion perpendicular to it such that *D*_{∥} > *D*_{⊥} as found in experiments.^{57} However, as one approaches the NSmA phase transition, diffusion along the mesogens’ axes is inhibited by the formation of smectic layers, such that *D*_{∥} ≈ *D*_{⊥} [see Figure 5(b)]. Deep in the smectic A phase at low temperature this effect is more pronounced and *D*_{∥} < *D*_{⊥}. Again this is in very good agreement with experimental observations.^{57}

As the density increases steadily with decreasing *T*, both diffusion coefficients become rather small. Because the smectic layers become more robust with decreasing *T*, *D*_{∥} nearly vanishes at the lower end of the temperature range considered in Figure 5(a) whereas *D*_{⊥} remains relatively substantial reflecting the two-dimensional fluidic character of the smectic layers.

At this stage we would like to emphasize that data shown in Figures 2 and 5 are fully consistent with each other and give a complete picture of the NSmA phase transition. The NSmA phase transition appears rather rounded in comparison with the corresponding IN phase transition in Figure 2. In the vicinity of the NSmA phase transition the formation of smectic layers is enhanced. The layers are not equidistant: vary in shape and size. Because of the finite size of the system the formation of irregular layers does not average out and results in a non-vanishing value of Λ. However, in the thermodynamic limit the effect vanishes. Therefore, the transition from *D*_{∥} > *D*_{⊥} to *D*_{∥} < *D*_{⊥} appears in our model system slightly before the NSmA phase transition in comparison with experiments where this transition occurs right after the NSmA phase transition.^{57}

Before closing the present section on equilibrium properties, we finally consider elastic properties of our model liquid crystal. To that end, we follow earlier work by Allen and Frenkel^{58,59} and Allen *et al.*^{60} and introduce the Fourier transform of the alignment tensor

where Eq. (4.1) has also been used. Let $e\u03021$, $e\u03022$, and $e\u03023$ be the axes of a coordinate system in which **Q** is diagonal and $n\u03020=0,0,1T$. For a given wave vector $k=k1,0,k3$, we introduce

where $Q\u02dc\alpha \beta $ (*α*, *β* = 1, 2, 3) are components of $Q\u02dc$ and *K*_{1}, *K*_{2}, and *K*_{3} are the Frank elastic constants. The above expressions are valid if the deviation of the director field $n\u0302$ from $n\u03020$ occurs on a length scale that is large compared to the correlation length so that we are focusing on the limit of small wave numbers *k*_{1}, *k*_{2}, and *k*_{3}. We refer the reader interested in a derivation of expressions presented in Eqs. (4.10a) and (4.10b) to the Appendix.

Plots in Figure 6 reveal the anticipated linear variation of *E*_{13} and *E*_{23} with $k12$ and $k32$ in the limit of sufficiently small wave numbers. However, as *T* is lowered towards the NSmA phase transition the simulation data exhibit a substantial scatter. This is because fluctuations of the Fourier components $Q\u02dc\alpha \beta $ become exceedingly small at lower *T*.

From the slopes of the curves plotted in Figure 6 we obtain the Frank elastic constants *K*_{1}, *K*_{2}, and *K*_{3}. The plots in Figure 7 reveal that the so-called one-constant approximation *K*_{1} = *K*_{2} = *K*_{3} frequently employed in theoretical studies of liquid-crystalline materials^{16} is approximately valid for sufficiently high temperatures *T* ≳ 0.82. Below this threshold and for decreasing *T*, all three elastic constants are different and hence the one-constant approximation is no longer satisfied. In particular, *K*_{2} and *K*_{3} presumably diverge right at the NSmA phase transition whereas *K*_{1} seems to remain finite all the way down to that transition. This is in good qualitative agreement with the experimental data of Madhusudana and Pratibha for the liquid crystal 8CB.^{61}

Moreover, the temperature dependence of the twist and bend elastic constants follows an exponential law, i.e., *K*_{2}, *K*_{3} ∝ (*T* − *T*_{NSmA})^{−ν} which was predicted by de Gennes^{62} theoretically and confirmed in three independent experiments.^{63–65} The exponent *ν* describing the critical behavior of the elastic constants is 0.50 in the framework of a mean field approach and 0.66 for experimental systems. Fitting the power law to our data [see Figure 7(b)], we obtain *ν* = 0.55 which is between the mean field model and an experimental system. The anomaly of *K*_{2} and *K*_{3} in the nematic phase is rationalized by the formation of cybotactic clusters. Cybotactic clusters reported by de Vries^{66} are smectic droplets which arise in the vicinity of the NSmA phase transition due to fluctuations in local density.^{62} The qualitative dependence of the Frank elastic constants on *T* as well as their critical behavior lends additional credibility to our model system. At this point it is noteworthy that the rotational viscosity, defined as the ratio of viscous torque and the resulting angular velocity, diverges in the vicinity of the NSmA phase transition as well (not shown here).^{67,68}

### B. Nonequilibrium steady states of a smectic A phase

Among the quantities that we wish to investigate by means of NEMD simulations is the velocity field generated at steady state by the external body force *F*_{e}. Using essentially macroscopic arguments de Gennes and Prost^{16} argue that in the direction of flow (i.e., in the *x*-direction for our setup)

where each prime denotes a derivation with respect to the *z* coordinate, $\eta \xaf$ is a measure of viscosity (in fact, as stated by de Gennes and Prost, it is “a certain average of the Leslie coefficients”^{16}), and *κ* is the inverse thickness of a boundary layer near the substrates in which the flow of the liquid crystal is dominated by friction. In essence, Eq. (4.11) is a statement about force balance in a steady-state nonequilibrium situation.

It should also be noted that the left hand side of Eq. (4.11) holds only if one assumes the phase under flow to be a homogeneous continuum.^{69} Whereas, strictly speaking, this is not exactly valid at the molecular level, we tacitly ignore the inhomogeneity of the smectic phase on account of the macroscopic character of the differential equation given in Eq. (4.11).

This equation can be easily solved by standard techniques and the general solution can be cast as $vx=\rho Fe/\eta \xaf\kappa 2+C1e\kappa z+C2e\u2212\kappa z$, where *C*_{1} and *C*_{2} are two constants of integration. Using the boundary condition $vzsz/2=vz\u2212sz/2=0$, it turns out that only the sum *C*_{1} + *C*_{2} = *C* matters in this problem where

and hence

solves the present boundary value problem.

Clearly, in this description of the velocity profile, *κ* introduces a characteristic length scale. If, on the one hand, *κ* is of a typical molecular magnitude, Eq. (4.13) describes plug flow. That plug flow is possible in a steady-state nonequilibrium situation and arises in various ordered liquid-crystalline mesophases has been noted by Helfrich a long time ago.^{19} He bases his analysis on essentially macroscopic arguments as we do here, too.

If, on the other hand, the friction dominated boundary layer becomes macroscopic in size (i.e., for *κ* → 0), the term in brackets in Eq. (4.13) can be recast as

Inserting this last expression into Eq. (4.13) we obtain

which is the characteristic parabola describing the velocity profile in Poiseuille flow. This expression has been derived earlier by Todd *et al.* under different but comparable assumptions.^{69–71} Clearly, *v*_{x} is largest at the midpoint between the substrates (*z* = 0) and vanishes directly at them (*z* = ± *s*_{z}/2) as it must. It is also a simple matter to verify that the far right side of Eq. (4.15) can be obtained directly from Eq. (4.11) by neglecting the term proportional to *κ*^{2}.

However, as we shall see shortly, the velocity profile obtained here is not exactly characteristic of Poiseuille flow. This is because the assumption of a diverging boundary-layer thickness (i.e., *κ* → 0) is only a rather rough approximation here. Therefore, the more detailed analysis of the velocity profiles below will always be based upon Eq. (4.13) rather than Eq. (4.15). However, we notice in passing that in the work of Stieger *et al.*,^{72} where a nematic rather than a smectic phase is considered, the far right hand side of Eq. (4.15) is valid and provides an excellent description of the velocity profile.

To demonstrate that the above macroscopic treatment can be applied at the microscopic level as well, we present in Figure 8 plots of *v*_{x} versus *z*. In all plots the flow direction is parallel to $e\u0302x$. As indicated by the cartoons in Figures 8(a) and 8(b), the mesogens are anchored at the substrates in a direction orthogonal and parallel to the flow direction, respectively. Irrespective of $n\u03020$ (i.e., the normal to the smectic layers, see also below), mesogens in these layers are always aligned with the plane of the solid substrates. Hence, the normal to the smectic A layers points in a direction orthogonal to the applied flow [see Figure 8(a)] or parallel with it as well [see Figure 8(b)].

Because of this orientation of the smectic layers, one observes Poiseuille-like flow under the conditions illustrated in Figure 8(a). Indeed, the parabolic profile predicted by Eq. (4.15) is obtained in the NEMD simulations except for the immediate vicinity of the solid substrates where the fit of Eq. (4.15) is not entirely capable of describing the NEMD data. Nonetheless, the overall quality of the fit is still surprisingly good indeed.

If instead one uses Eq. (4.13), plots in Figure 8(a) reveal that this equation describes the NEMD data almost perfectly across the entire space between the solid substrates. This is also true for data shown in Figure 8(b) where the different orientation of the smectic layers gives rise to plug instead of Poiseuille-like flow. The characteristic feature of plug flow, namely *v*_{x} ≃ const, holds quite well for a large part of the LC around its midsection.

However, the thickness of the boundary layer *κ*^{−1} diverges by no means but remains even smaller than *s*_{z} which puts an upper limit to *κ*^{−1}. This is illustrated by the fact that the fit of Eq. (4.13) to the NEMD data yields *κ*^{−1} ≃ 5.56 in the case of data shown in Figure 8(a); *κ*^{−1} ≃ 2.44 is obtained from the fit in Figure 8(b).

Hence, the thickness of the boundary layer remains on a molecular scale and, therefore, the far right hand side of Eq. (4.15) constitutes a rather rough approximation to the NEMD generated velocity profiles. In other words, the flow observed in the case in which the plane of the smectic layers is orthogonal to the direction of flow is only Poiseuille-like. Nonetheless, it seems worthwhile emphasizing that a macroscopic picture of the flow underlying Eq. (4.11) provides an excellent description under the conditions of the present study despite our otherwise molecular approach.

We note in passing that if the mesogens are anchored homeotropically such that the smectic layers form in the *x*–*y* plane and are stacked along the *z*-axis, the system cannot reach a steady nonequilibrium state if flow is initiated in the *x*- (or *y*-) direction. This is because mesogens in neighboring layers do not experience any frictional forces acting between them. In fact, the plot in Figure 1 illustrates this quite nicely because the equipotential surface is purely repulsive for end-end configurations of a mesogenic pair. Consequently, entire smectic layers flow more or less independently past each other as individual entities such that *v*_{x} changes discontinuously from one layer to the next. The velocity profile rises towards the midpoint of the liquid crystal. This is because as one moves away from a solid substrate, mesogens also experience less friction from that substrate.

We are, however, in a position to extract more detailed information from flow profiles such as the ones plotted in Figures 8(a) and 8(b). We begin this analysis by presenting in Figure 9 plots of the viscosity $\eta \xaf$ as a function of the Péclet number^{15}

where

is the mean velocity of the flow and $D\xaf\u226112(D\u2225+D\u22a5)$. Because of its definition in Eq. (4.16), the Péclet number expresses the relative contribution of advective to diffusive motion of the mesogens. Regardless of whether one considers Poiseuille-like or plug flow, $\eta \xaf$ increases monotonically with ℘. As one can see from Figure 9, $\eta \xaf$ for a Poiseuille-like flow situation exceeds that under plug-flow conditions.

Generally speaking, an increase in viscosity reflects an increase in friction in the system. That this is lower in plug flow is reflected by plots in Figure 8(b) which indicate that a relatively large portion of the smectic phase is moving more or less at the same velocity. Intuitively, one then expects friction experienced by mesogens located in this region to be lower than in Poiseuille-like flow where *v*_{x} varies all across the smectic phase as revealed by Figure 8(a). Hence, plots in Figure 9 are consistent with data presented in both parts of Figure 8.

The difference in viscosity can be ascribed to changes in the structure of the liquid crystal under flow conditions. This can be rationalized through plots of the nematic and smectic order parameters *S* and Λ in both parts of Figure 10. One can see that in a Poiseuille-like flow both order parameters appear to be independent of ℘ over the entire range considered [see Figure 10(a)] whereas in plug flow [see Figure 10(b)], *S* declines steadily until at ℘ = 1550 it reaches only about 60% of its value in Poiseuille-like flow. At the same time, Λ remains constant and is independent of the specific flow conditions as a comparison of plots in Figures 10(a) and 10(b) clearly indicates.

The decline of *S*, which is a *globally* defined quantity, implies that the director field $n\u0302$ is perturbed *locally*. We compute $n\u0302$ by solving an eigenvalue equation similar to Eq. (4.4) where, however, **Q**_{V} is replaced by its local counterpart **Q** from Eq. (4.1). We present plots of $n\u0302$ in Figures 11(a)–11(c).

In principle, three geometrically distinct deformations of $n\u0302$ would be possible in an LC composed of uniaxial mesogens, namely splay, twist, and bend. Associated with each one of these is a contribution to the Frank elastic free energy with the respective coupling constants *K*_{1}, *K*_{2}, and *K*_{3}. Because of our results already presented in Figure 7, we see that upon approaching the NSmA phase transition *K*_{2} and *K*_{3} increase (and may even diverge^{61}) whereas *K*_{1} remains smallest (and finite). Thus, for the LC to optimize its elastic free energy, it seems best to avoid twist and bend deformations altogether.^{16,73}

Assuming now splay to be the only “surviving” deformation of $n\u0302$ and starting from equidistant smectic layers in the LC at rest, immediately leads one to conclude that this equidistance will be preserved in a nonequilibrium steady-state situation. This is indeed the case as Figure 11(d) and the plot of Λ in Figure 10(b) reveal. Moreover, the plot also shows that the periodic sequence of planar smectic layers in an LC at rest is transformed into concave planes with respect to the direction of flow in nonequilibrium steady-state situations, where the curvature of each plane is the same at the same *x* and increases with ℘.

Because of this curved structure of smectic layers under flow, we realize that around the midsection of the LC centered on *z* = 0 mesogens may still align with the direction of flow which coincides with the direction of wall anchoring. Hence, for mesogens located in a small volume centered on *z* = 0, $cos\alpha \u2261u\u0302\omega \u22c5e\u0302z=0$ ($\alpha =\pi 2$) on average. Likewise, $\alpha =\pi 2$ for mesogens in the immediate vicinity of the walls on account of anchoring [see Eqs. (2.8) and (2.9)]. Starting now at *z* = 0 and moving along the positive *z*-axis, the curved structure of the smectic layers requires the deflection angle *α* to decrease from its initial value of $\pi 2$ and then to increase again as the upper wall is approached where $\alpha =\pi 2$ is once again attained. The same is expected as one approaches the lower wall starting at the system’s midpoint *z* = 0 only that in this case *α* *increases* at first.

From these considerations, it becomes apparent that under flow one anticipates two parallel bands to exist along the flow direction in which the director field is deformed the way just described. This is indeed the case as plots in Figures 11(a)–11(c) clearly show.

To quantify the deformation of $n\u0302$, we introduce the elastic Frank free-energy density associated with only splay deformations of $n\u0302$ via

We note in passing that in general an additional term would arise in Eq. (4.18) affiliated with the spacing between smectic layers. However, if this distance is constant, as we have assumed from the outset, this additional contribution to the elastic free energy vanishes and can therefore safely be ignored in Eq. (4.18).^{73} Plots in Figures 11(a)–11(c) show that as the flow is cranked up (i.e., with increasing ℘) the elastic Frank free-energy density associated with the deformation bands of $n\u0302$ increases which implies a stronger deformation of $n\u0302$ as the strength of the flow increases.

From a slightly different perspective, one can also relate the variation of the deflection angle *α* to the shear rate $\gamma \u0307\u2261\u2202vx/\u2202z$. As one can see from Figure 12, at sufficiently large $\gamma \u0307$, $\alpha \u2192\pi 2$. These shear rates are obtained in regions sufficiently close to the walls where the competition between anchoring-induced alignment of the mesogens and the flow-induced distortion of $n\u0302$ is largest. With decreasing $\gamma \u0307$, the deflection angle decreases (increases) until a minimum (maximum) is reached. These extrema occur near the middle of the distortion bands visible in particular in Figures 11(b) and 11(c). A further decrease of $\gamma \u0307$ then causes *α* to quickly approach $\pi 2$ corresponding to a perfect flow alignment of the mesogens at the midsection of the LC where, consequently, $\gamma \u0307$ vanishes. Notice also that the variation of *α* with $\gamma \u0307$ in the small shear-rate regime is more pronounced at higher Péclet numbers. This is consistent with the plots in Figures 11(b) and 11(c) which show that the distortion of the director field is stronger the larger ℘ is.

The apparent insensitivity of Λ to the flow illustrated in Figure 10(b) can now be understood from the plots in Figure 11. Notice that the local director field is unaffected by the flow around the midsection of the liquid crystal. Moreover, it is symmetric with respect to this midsection such that its bending in the upper part is compensated by that in its lower part. As a consequence, $n\u03020$ stays fixed regardless of the strength of the flow. In addition, the periodicity of the local density is preserved everywhere. Both features, a fixed $n\u03020$ and the preserved periodicity of the smectic layers, cause Λ to remain unaltered as ℘ increases.

Nonetheless, all three plots in Figure 11 clearly show that a lot of nematic order is still preserved *locally* which is the reason why *S* in Figure 10(b) remains relatively large even at ℘ = 1550. Moreover, it is worth pointing out that the local nematic order parameter remains homogeneous for all the cases considered in Figure 11 (and is consequently not shown).

## V. DISCUSSION AND CONCLUSIONS

In this work we present extensive equilibrium and nonequilibrium computer simulations of a simple model liquid crystal that allows for the formation of nematic as well as smectic A phases. In this model, mesogens have an aspect ratio exceeding that of a previously considered model version.^{36} Even though this aspect ratio is still much smaller than for mesogens in the more widely used Gay-Berne model, it turns out that the elongated shape of the mesogens is a necessary prerequisite for the formation of a smectic A phase. Under thermodynamic conditions chosen in this work we observe a stable nematic phase over a temperature range 0.87 ≲ *T* ≲ 0.78; above and below this range isotropic and smectic A phases are stable, respectively.

By computing mean square displacements in directions parallel and perpendicular to the mesogens’ long axes, we investigate the dynamics of our model liquid crystal at thermodynamic equilibrium. In the limit of sufficiently long times the mesogens exhibit ordinary diffusive behavior. However, an intermediate time scale exists where self-diffusion in the direction normal to the plane of the smectic layers is subdiffusive. This is a reflection of the one-dimensional crystal-like structure of a smectic A phase. Once a mesogen reaches the boundary of its original layer, diffusion is hindered on account of the presence of the next layer. It then takes a while until mesogens located in that next layer make room such that a new mesogen can enter that layer. Eventually, at sufficiently long times this will happen because of the two-dimensional fluidic structure of each smectic layer.

As one might have anticipated from the above observations, self-diffusion is anisotropic in the nematic and smectic A phases. This anisotropy is manifested through the relation between *D*_{∥} and *D*_{⊥}. In the nematic phase *D*_{∥} > *D*_{⊥}, whereas *D*_{∥} < *D*_{⊥} in the smectic phase. Similar behavior is found in real systems using NMR diffusometry.^{57}

A feature of our model that might seem to indicate that it is suitable for describing cyanobiphenyl-based liquid crystalline materials is indicated by its elastic behavior. We compute the elastic constants for splay, twist, and bend deformations of the director field at thermodynamic equilibrium from fluctuations in the Fourier components of the local alignment tensor.

Even though these data exhibit a substantial amount of scatter (especially at low wave numbers) as one approaches the NSmA phase transition by lowering the temperature, they still permit to extract the elastic constants reliably. If the latter are plotted as functions of *T* in the nematic phase, we observe that the elastic constant for splay deformations of the director field apparently remains finite all the way to the NSmA phase transition; the other two elastic constants seem to diverge instead, given our numerical resolution. This behavior has also been seen experimentally for octylcyanobiphenyl.^{61}

The simulations of static and dynamic properties of our model at thermodynamic equilibrium are complemented by nonequilibrium steady-state simulations. Here the smectic A phase is particularly interesting on account of its one-dimensional solidlike character. The broken symmetry in the direction normal to the plane of the smectic layers permits one to initiate flow either parallel with the plane of the smectic layers or perpendicular to this plane. In the former case one observes a Poiseuille-like velocity profile whereas in the latter situation plug flow is observed.

By using macroscopic arguments borrowed from de Gennes and Prost’s book,^{16} we see that for a smectic A phase there are small but significant deviations from perfect Poiseuille flow very much akin to the case where the liquid crystal is in its nematic phase.^{72} If flow is initiated in a direction perpendicular to the plane of the smectic layers, a large portion of the mesogens travel at the same speed which is a characteristic of plug flow. This was shown many years ago—again from a macroscopic perspective—by Helfrich^{19} who pointed out the possibility of plug flow in cholesteric and smectic liquid crystals.

Comparing Poiseuille-like flow with plug flow, it turns out that the viscosity is lower in the latter case. We believe that this is because in the portion of the liquid crystal in which mesogens are travelling at the same speed, they experience less friction from neighboring portions of the fluid compared with Poiseuille-like flow where no such portion of the liquid crystal exists.

Last but not least, a somewhat surprising observation concerns the structure of the liquid crystal in plug flow. Over the range of Péclet numbers accommodated in this work, the layer structure characteristic of a smectic A phase at thermodynamic equilibrium is preserved overall but bent locally. In fact, the smectic layers form a periodic sequence of concave surfaces with respect to the direction of flow.

## Acknowledgments

We are grateful for financial support from the Deutsche Forschungsgemeinschaft via the International Graduate Research Training Group 1524. One of us (R.G.) thanks the *Studienstiftung des Deutschen Volkes* for financial support. Another one (C.K.H.) gratefully acknowledges support from the US National Science Foundation under Grant No. OISE 1065466 and the Research Triangle MRSEC under Grant No. DMR-1121107.

### APPENDIX: ELASTIC CONSTANTS

For the convenience of the reader we derive the expressions in Eqs. (4.10a) and (4.10b) which have been used to measure the director’s fluctuations in order to determine Frank elastic constants. Therefore we follow Oswald and Pieranski^{74} and cast the Frank free energy according to

where the terms on the right-hand side represent contributions of splay, twist, and bend deformations, respectively. Perturbations of the global director, which is considered to be $n\u03020=(0,0,1)T$, occur perpendicular to it. Hence, the local director may be cast as

where $\delta n\u0302=(\delta n1,\delta n2,0)T$ is a perturbation. For the sake of convenience we focus on the splay deformation of the free energy. Because $n\u03020$ is pointing in the *z*-direction we may rewrite the splay contribution to the frank free energy according to

where ∂_{α}*n*_{β} ≡ ∂*n*_{β}/∂*α* and *α*, *β* = 1, 2, 3. The Fourier transform (indicated by the tilde) of the components of the local director $n\u0302(r)$ is given by

where ** k** is the wave vector. Consequently, the inverse transformation is obtained via

which we substitute in Eq. (A3) and obtain

This expression for $Fsplay$ can be easily simplified by the application of the well known relations

where *δ* is the Dirac *δ* distribution and

such that

Furthermore, using $n\u02dc\alpha \u2009(\u2212k)=n\u02dc\alpha \u2217\u2009(k)$ and $n\u02dc\alpha \u2009(k)n\u02dc\alpha \u2217\u2009(k)=|n\u02dc\alpha (k)|2$ (the asterisk denotes the complex conjugate), and replacing the integral by a sum over the wave vectors, we cast the final result as

In the same manner we obtain

and

for the free energy of the twist and bend deformations, respectively. We define a coordinate system **1**, **2**, **3**, where the global director is pointing in the *z*-direction and the wave vector is chosen to be in the 1-3 plane, such that ** k** = (

*k*

_{1}, 0,

*k*

_{3})

^{T}. Furthermore, we take the unit vector $e\u03021$ to be parallel to the plane spanned by

**and $n\u03020$ and $e\u03022$ perpendicular to it. Thus, the wave vector $k=k1e\u03021+k3n\u03020$ is used to rewrite and summarize Eqs. (A10)–(A12) to obtain**

*k*It is helpful to remember that $n\u02dc\alpha (k)$ is the Fourier transform of the components of $n\u0302(r)$ which describes local fluctuations of the director field. The fluctuations of the director field occur due to molecular rotation. Because all three rotational degrees of freedom are decoupled, we use the equipartition theorem to separate the director field fluctuations according to

by analogy with Allen and Frenkel^{58,59} who used the tensorial version for their derivation. In order to rewrite the expressions presented in Eqs. (A14) and (A15) as their tensorial analogs, we use the macroscopic version of the alignment tensor [see Eq. (4.1)]

and remember that $n\u0302(r)=(\delta n1,\delta n2,1)T$ [see Eq. (A2)]. Substituting $n\u0302(r)$ in Eq. (A16) results in $Q13(r)=32S\delta n1$ and $Q23(r)=32S\delta n2$, where the computation of *δn*_{1} and *δn*_{2} is performed in Fourier space via Eqs. (A14) and (A15). Using this proportionality one immediately can rewrite Eqs. (A14) and (A15) in Eqs. (4.10a) and (4.10b), respectively.

## REFERENCES

Gerolamo Cardano (1501–1576), Italian mathematician, physician, inventor, astrologer, and gambler.