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 K1, K2, and K3 for the respective splay, twist, and bend deformations of the director field n̂ 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.

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 microlasers1–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 analogues7,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 grown15 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 Pershan20 and later considered in simplified form by de Gennes.21 More recently, Bennett and Hess22 studied the Miesowicz, Helfrich, and Leslie viscosities in proximity of the nematic to smectic phase transition. Walker and Stewart23 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 equilibrium31,32 and nonequilibrium33,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.

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

ΦmmR,Ω=i=1N1j=i+1Nφmmrij,ωi,ωj,
(2.1)

where rij = rirj is the distance vector connecting the centers of mass of mesogens i and j located at ri and rj, respectively. In addition to rij, 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,,rN and Ω=ω1,,ωN to indicate that Φmm depends on the configuration of the N mesogens. On account of the uniaxial symmetry of the mesogens, ω=θ,ϕ 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

φmmrij,ωi,ωj=φisorij+φanisrij,ωi,ωj,
(2.2)

where rij=rij. More specifically, we take

φisorij=εa1σrij10a2expλrijrij=φreprij+φattrij,
(2.3)

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

a1=1101+λσ9λσ,
(2.4a)
a2=σexpλσ9λσ
(2.4b)

are introduced to guarantee that the minimum of φiso remains at rijmin=σ and that φiso(rijmin)=ε/10 regardless of λ.

For the anisotropic interactions, we follow Giura and Schoen38 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 ωiωi=ωi or ωjωj=ωj as well as its invariance upon replacing rij by −rij. One may then cast φanis as

φanisrij,ωi,ωj=φattrijΨr̂ij,ωi,ωj,
(2.5)

where r̂ij=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

Ψr̂ij,ωi,ωj=ε1P2ûωiûωj+ε2P2ûωir̂ij+P2ûωjr̂ij
(2.6)

where P2x=123x21 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 û(ωi)û(ωj)=±1 (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., φiso(rijmin)=ε/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 ε1<ε2, φanis is attractive on the one hand as long as û(ωi)û(ωj)=±1and in addition û(ωi)r̂ij=û(ωj)r̂ij=0. If, on the other hand, û(ωi)û(ωj)=û(ωi)r̂ij=û(ωj)r̂ij=±1, φ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.

FIG. 1.

Contour plot of φmm [see Eqs. (2.2)–(2.6)] for a pair of mesogens whose centers of mass are located in the x-y plane such that r12 = (x12, y12, 0)T (T denotes the transpose). Both mesogens have the orientation û(ω1)=û(ω2)=(1,0,0)T. The color bar indicates the value of φmm and the black lines are equipotential lines along which φmm assumes the values given in the figure.

FIG. 1.

Contour plot of φmm [see Eqs. (2.2)–(2.6)] for a pair of mesogens whose centers of mass are located in the x-y plane such that r12 = (x12, y12, 0)T (T denotes the transpose). Both mesogens have the orientation û(ω1)=û(ω2)=(1,0,0)T. The color bar indicates the value of φmm and the black lines are equipotential lines along which φmm assumes the values given in the figure.

Close modal

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 chemicals40,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 sz along the z-axis. The contribution of the mesogen-substrate interaction to the total configurational potential energy can then be cast as

ΦmsR,Ω;S=i=1Nj=12Nsφmsrij,ωi
(2.7)

assuming each substrate to be composed of a monolayer of Ns 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 /σ=43 corresponding to an areal density of the substrate atoms of ρsσ2 = 2/ℓ2. Taking the configuration of substrate atoms S = (s1, …, s2Ns) to be fixed, rij=risj 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),

φmsrij,ωi=εσrij12σrij6gαωi,
(2.8)

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

gαωi=ûωiêα2,
(2.9)

where êα (α = 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 û is parallel to the versor; on the contrary, if û 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 2NNs 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̂0 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.

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 VV′, RR=RV/V3 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 rc = 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 rN = 3.8σ. This reduces the computational cost associated with the search for mesogens that are within the radius rc 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 × 103 such cycles for equilibration followed by another 5 × 104 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 ri 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 thermostat48–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 rc and rN are used. We equilibrate the system for 2.0 × 105 to 3.0 × 105 time steps followed by another 106 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 × 104 mesogens placed in a rectangular nanofluidic channel of volume V = sxsysz, 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.

FIG. 2.

Plot of the nematic S (red square) and smectic order parameter Λ (blue circle) as functions of temperature T. (Black circle) demarcates the IN phase transition (see text). Lines are fits intended to guide the eye.

FIG. 2.

Plot of the nematic S (red square) and smectic order parameter Λ (blue circle) as functions of temperature T. (Black circle) demarcates the IN phase transition (see text). Lines are fits intended to guide the eye.

Close modal

In the NEMD simulations the equations of motion are integrated numerically using again the velocity Verlet algorithm proposed by Ilnytskyi and Wilson47 with a time step of δt = 10−3. A steady-state flow is initiated and maintained by applying a body force Fe=Feêx 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 ε/kB, where kB is Boltzmann’s constant. Other quantities such as time, pressure, or density can be expressed in terms of combinations of the basic ones.46 

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 Frenkel51 and introduce the alignment tensor

Qr=12ρri=1N3ûωiûωi1δrri
(4.1)

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

ρr=i=1Nδrri
(4.2)

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 indicates an ensemble average in the isothermal isobaric ensemble.

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

QV=12Ni=1N3ûωiûωi1
(4.3)

which satisfies the eigenvalue equation

QVn̂k=λkn̂k
(4.4)

and where n̂k is the kth eigenvector corresponding to eigenvalue λk.

Because QV 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’s52 formula.53 In an infinitely large system the discriminant D3 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, D3 < 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=QV gives the ensemble averaged eigenvalues λk, where in general λk=λk. Notice, however, that on account of the head-tail symmetry of the mesogens, n̂kn̂k. Following Eppenga and Frenkel we then define the nematic order parameter via S=maxk=13λk and take the associated eigenvector as the nematic director n̂0.

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

Λd=1Nj=1Nexp2πirjn̂0d.
(4.5)

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̂0. 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,Ω 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 Schoen37 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, rjn̂0/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, rjn̂0/d is not an integer but a real number in the absence of smectic layers (for example, in the nematic phase where n̂0 is still well defined). If in this case the limit N → ∞ is considered, j=101dxexp2π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 TIN 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

g20λ02λ02.
(4.6)

Greschek and Schoen37 found and subsequently rationalized that for a sufficiently weak discontinuous phase transition g201 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 g20N.37 Thus, the larger N the larger is g201. 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 TIN ≃ 0.87 of the IN phase transition.

FIG. 3.

Plots of the second-order Binder cumulant g201 as a function of temperature T and for various system sizes: N = 500 (red square), N = 1000 (blue circle), N = 2500 (black triangle), N = 5000 (downward magenta triangle). The horizontal line at the intersection of the curves marks the temperature TIN ≃ 0.87 at the IN phase transition. Inset is an enlargement in the vicinity of TIN and lines are fits intended to guide the eye.

FIG. 3.

Plots of the second-order Binder cumulant g201 as a function of temperature T and for various system sizes: N = 500 (red square), N = 1000 (blue circle), N = 2500 (black triangle), N = 5000 (downward magenta triangle). The horizontal line at the intersection of the curves marks the temperature TIN ≃ 0.87 at the IN phase transition. Inset is an enlargement in the vicinity of TIN and lines are fits intended to guide the eye.

Close modal

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 g201 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, TIN from the cumulant analysis agrees quite well with the estimated inflection point of the curve ST plotted in Figure 2. At TIN, SIN ≃ 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, TNSmA ≃ 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=riûωiûωi and ri=riri. With these quantities we introduce the mean-square displacements (MSDs),

Δr,2t=1Ni=1Nri,t0+tri,t02t0,
(4.7)

where t0 indicates an average over a set of different initial times t0 due to stationarity56 of the self-diffusion process. With Eq. (4.7) we define

D=limtΔr2t2t,
(4.8a)
D=limtΔr2t4t
(4.8b)

as quantitative measures of self-diffusion at thermodynamic equilibrium.

Plots of Δr2 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, Δr2t2. 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.

FIG. 4.

(a) Plots of the mean square displacements Δr2t as functions of time t on a double-logarithmic scale. (b) As (a), but for Δr2t. In both parts of the figure (red square) T = 0.90 (isotropic), (blue circle) T = 0.84 (nematic), (black triangle) T = 0.70 (smectic A) (see also Figure 2). The magenta solid line indicates a slope proportional to t.

FIG. 4.

(a) Plots of the mean square displacements Δr2t as functions of time t on a double-logarithmic scale. (b) As (a), but for Δr2t. In both parts of the figure (red square) T = 0.90 (isotropic), (blue circle) T = 0.84 (nematic), (black triangle) T = 0.70 (smectic A) (see also Figure 2). The magenta solid line indicates a slope proportional to t.

Close modal

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 Δr2t 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., Δr2tν, ν < 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̂0 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, DD which makes sense because no spatial direction is distinguished here.

FIG. 5.

(a) Semi-logarithmic plots of diffusion coefficients D (red square) and D (blue circle) as functions of temperature T. (b) Plots of the ratio D/D (downward magenta triangle) as a function of T. Vertical lines (black solid line and black dashed line) mark the IN and NSmA phase transitions, respectively (see Figure 2).

FIG. 5.

(a) Semi-logarithmic plots of diffusion coefficients D (red square) and D (blue circle) as functions of temperature T. (b) Plots of the ratio D/D (downward magenta triangle) as a function of T. Vertical lines (black solid line and black dashed line) mark the IN and NSmA phase transitions, respectively (see Figure 2).

Close modal

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 DD [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 Frenkel58,59 and Allen et al.60 and introduce the Fourier transform of the alignment tensor

Q˜k=drQrexpikr,
(4.9)

where Eq. (4.1) has also been used. Let ê1, ê2, and ê3 be the axes of a coordinate system in which Q is diagonal and n̂0=0,0,1T. For a given wave vector k=k1,0,k3, we introduce

E13k1,k394S2VkBT|Q˜13|2=k1,k30K1k12+K3k32,
(4.10a)
E23k1,k394S2VkBT|Q˜23|2=k1,k30K2k12+K3k32,
(4.10b)

where Q˜αβ (α, β = 1, 2, 3) are components of Q˜ and K1, K2, and K3 are the Frank elastic constants. The above expressions are valid if the deviation of the director field n̂ from n̂0 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 k1, k2, and k3. 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 E13 and E23 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˜αβ become exceedingly small at lower T.

FIG. 6.

Plots of (a) E13k1, (b) E23k1, and (c) E13k3 as functions of the wavenumbers k12 or k32 and for temperatures T = 0.86 (red square), T = 0.84 (blue circle), T = 0.82 (black triangle), and T = 0.80 (downward magenta triangle), all of which are above the NSmA phase transition. Lines represent linear fits to the discrete simulation data.

FIG. 6.

Plots of (a) E13k1, (b) E23k1, and (c) E13k3 as functions of the wavenumbers k12 or k32 and for temperatures T = 0.86 (red square), T = 0.84 (blue circle), T = 0.82 (black triangle), and T = 0.80 (downward magenta triangle), all of which are above the NSmA phase transition. Lines represent linear fits to the discrete simulation data.

Close modal

From the slopes of the curves plotted in Figure 6 we obtain the Frank elastic constants K1, K2, and K3. The plots in Figure 7 reveal that the so-called one-constant approximation K1 = K2 = K3 frequently employed in theoretical studies of liquid-crystalline materials16 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, K2 and K3 presumably diverge right at the NSmA phase transition whereas K1 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 

FIG. 7.

(a) Plots of the Frank elastic constants K1 (red square), K2 (blue circle), and K3 (black triangle) as functions of temperature T. Lines are fits to guide the eye. (b) Same as in (a) but here blue solid line and black solid line are fits of K(T) = a(TTNSmA)ν to the twist and bend elastic constants, respectively (see text).

FIG. 7.

(a) Plots of the Frank elastic constants K1 (red square), K2 (blue circle), and K3 (black triangle) as functions of temperature T. Lines are fits to guide the eye. (b) Same as in (a) but here blue solid line and black solid line are fits of K(T) = a(TTNSmA)ν to the twist and bend elastic constants, respectively (see text).

Close modal

Moreover, the temperature dependence of the twist and bend elastic constants follows an exponential law, i.e., K2, K3 ∝ (TTNSmA)ν which was predicted by de Gennes62 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 K2 and K3 in the nematic phase is rationalized by the formation of cybotactic clusters. Cybotactic clusters reported by de Vries66 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

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 Fe. Using essentially macroscopic arguments de Gennes and Prost16 argue that in the direction of flow (i.e., in the x-direction for our setup)

ρFeη¯=κ2vxzvxz,
(4.11)

where each prime denotes a derivation with respect to the z coordinate, η¯ 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=ρFe/η¯κ2+C1eκz+C2eκz, where C1 and C2 are two constants of integration. Using the boundary condition vzsz/2=vzsz/2=0, it turns out that only the sum C1 + C2 = C matters in this problem where

C=ρFeη¯κ21coshκsz/2
(4.12)

and hence

vxz=ρFeη¯κ21coshκzcoshκsz/2
(4.13)

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

1coshκsz/2cosh(κsz/2)coshκz=1cosh(κsz/2)1+12!κ2sz24112!κ2z2+O(κ4)κ22cosh(κsz/2)sz24z2.
(4.14)

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

vxz=ρFe2η¯cosh(κsz/2)sz24z2κ0ρFe2η¯sz24z2
(4.15)

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, vx is largest at the midpoint between the substrates (z = 0) and vanishes directly at them (z = ± sz/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 vx versus z. In all plots the flow direction is parallel to êx. 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̂0 (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)].

FIG. 8.

Plots of the velocity component vx as a function of position z between the solid substrates. Flow is initiated parallel to êx; (a) smectic layers perpendicular to the direction of flow, (red circle) NEMD results, (red solid line) fit with Eq. (4.13), (red dashed line), fit with Eq. (4.15); (b) as (a), but for smectic layers oriented in the direction of flow where the blue squares represent the NEMD data and the full line is a fit of Eq. (4.13) to these data. Insets are cartoons illustrating the orientation of smectic layers with respect to the direction of flow and the solid walls (grey shaded rectangular areas). Data are obtained for Fe = 0.05.

FIG. 8.

Plots of the velocity component vx as a function of position z between the solid substrates. Flow is initiated parallel to êx; (a) smectic layers perpendicular to the direction of flow, (red circle) NEMD results, (red solid line) fit with Eq. (4.13), (red dashed line), fit with Eq. (4.15); (b) as (a), but for smectic layers oriented in the direction of flow where the blue squares represent the NEMD data and the full line is a fit of Eq. (4.13) to these data. Insets are cartoons illustrating the orientation of smectic layers with respect to the direction of flow and the solid walls (grey shaded rectangular areas). Data are obtained for Fe = 0.05.

Close modal

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 vx ≃ 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 sz 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 xy 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 vx 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 η¯ as a function of the Péclet number15 

=szv¯xD¯,
(4.16)

where

v¯x=1szsz/2sz/2dzvxz
(4.17)

is the mean velocity of the flow and D¯12(D+D). 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, η¯ increases monotonically with ℘. As one can see from Figure 9, η¯ for a Poiseuille-like flow situation exceeds that under plug-flow conditions.

FIG. 9.

Plots of the viscosity η¯ as a function of the Péclet number ℘; (red square) Poiseuille flow [cf., Figure 8(a)]; (blue circle) plug flow [cf., Figure 8(b)].

FIG. 9.

Plots of the viscosity η¯ as a function of the Péclet number ℘; (red square) Poiseuille flow [cf., Figure 8(a)]; (blue circle) plug flow [cf., Figure 8(b)].

Close modal

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 vx 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.

FIG. 10.

Plots of the nematic and smectic order parameters S (red square) and Λ (blue circle) as functions of the Péclet number ℘; (a) Poiseuille-like and (b) plug flow.

FIG. 10.

Plots of the nematic and smectic order parameters S (red square) and Λ (blue circle) as functions of the Péclet number ℘; (a) Poiseuille-like and (b) plug flow.

Close modal

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

FIG. 11.

Maps of the director field n̂r represented by short dashes in the xz plane and for different strengths of the external body force Fe. Gray bars at the top and bottom of each frame represent the solid substrates and the arrows indicate the direction of flow; (a) ℘ = 127 (Fe = 0.01), (b) ℘ = 514 (Fe = 0.03), (c) ℘ = 1086 (Fe = 0.05) [cf., Figure 10(b)], (d) “snapshot” of a configuration from a NEMD simulation (Fe = 0.05). In part (d) mesogens are colored in blue if they are aligned with the direction of the flow. Thick orange lines mark the position of three bent smectic layers. In all parts of the figure gx has been used in Eq. (2.8) [see also Eq. (2.9)]. The attached color bar gives the local value of the Frank free energy density fd/K1 (see text).

FIG. 11.

Maps of the director field n̂r represented by short dashes in the xz plane and for different strengths of the external body force Fe. Gray bars at the top and bottom of each frame represent the solid substrates and the arrows indicate the direction of flow; (a) ℘ = 127 (Fe = 0.01), (b) ℘ = 514 (Fe = 0.03), (c) ℘ = 1086 (Fe = 0.05) [cf., Figure 10(b)], (d) “snapshot” of a configuration from a NEMD simulation (Fe = 0.05). In part (d) mesogens are colored in blue if they are aligned with the direction of the flow. Thick orange lines mark the position of three bent smectic layers. In all parts of the figure gx has been used in Eq. (2.8) [see also Eq. (2.9)]. The attached color bar gives the local value of the Frank free energy density fd/K1 (see text).

Close modal

In principle, three geometrically distinct deformations of n̂ 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 K1, K2, and K3. Because of our results already presented in Figure 7, we see that upon approaching the NSmA phase transition K2 and K3 increase (and may even diverge61) whereas K1 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̂ 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αûωêz=0 (α=π2) on average. Likewise, α=π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 π2 and then to increase again as the upper wall is approached where α=π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̂, we introduce the elastic Frank free-energy density associated with only splay deformations of n̂ via

fdrK1=12n̂r2.
(4.18)

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̂ increases which implies a stronger deformation of n̂ 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 γ̇vx/z. As one can see from Figure 12, at sufficiently large γ̇, απ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̂ is largest. With decreasing γ̇, 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 γ̇ then causes α to quickly approach π2 corresponding to a perfect flow alignment of the mesogens at the midsection of the LC where, consequently, γ̇ vanishes. Notice also that the variation of α with γ̇ 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.

FIG. 12.

Dependence of the deflection angle α (see text) on the shear rate, γ̇ (see text, cf., plots in Figure 11); (red square) ℘ = 514, (blue circle) ℘ = 1086. Inset is an enhancement for small shear rates.

FIG. 12.

Dependence of the deflection angle α (see text) on the shear rate, γ̇ (see text, cf., plots in Figure 11); (red square) ℘ = 514, (blue circle) ℘ = 1086. Inset is an enhancement for small shear rates.

Close modal

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̂0 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̂0 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).

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 Helfrich19 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.

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.

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 Pieranski74 and cast the Frank free energy according to

Fel=12VdrK1[n̂(r)]2+K2[n̂(r)×n̂(r)]2+K3[n̂(r)××n̂(r)]2,
(A1)

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̂0=(0,0,1)T, occur perpendicular to it. Hence, the local director may be cast as

n̂(r)=n̂0+δn̂,
(A2)

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

Fsplay=12K1Vdr[1n1(r)+2n2(r)]2,
(A3)

where ∂αnβ ≡ ∂nβ/∂α and α,  β = 1, 2, 3. The Fourier transform (indicated by the tilde) of the components of the local director n̂(r) is given by

n˜α(k)=drnα(r)exp(ikr),
(A4)

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

nα(r)=1(2π)3dkn˜α(k)exp(ikr)
(A5)

which we substitute in Eq. (A3) and obtain

Fsplay=12(2π)6K1drdkdk[n˜1(k)k1n˜1(k)k1+2n˜1(k)k1n˜2(k)k2n˜2(k)k2n˜2(k)k2]exp[i(k+k)r].
(A6)

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

drexp[i(k+k)r]=(2π)3δ(k+k),
(A7)

where δ is the Dirac δ distribution and

dkf(k)δ(k+k)=f(k)
(A8)

such that

Fsplay=12(2π)3K1dk[n˜1(k)n˜1(k)k12+2n˜1(k)k1n˜2(k)k2+n˜2(k)n˜2(k)k22].
(A9)

Furthermore, using n˜α(k)=n˜α(k) and n˜α(k)n˜α(k)=|n˜α(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

Fsplay=12VK1k[|n˜1(k)|k1+|n˜2(k)|k2].
(A10)

In the same manner we obtain

Ftwist=12VK2k[|n˜1(k)|k2|n˜2(k)|k1]
(A11)

and

Fbend=12VK3k[|n˜1(k)|2|n˜2(k)|2]k3
(A12)

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 = (k1, 0, k3)T. Furthermore, we take the unit vector ê1 to be parallel to the plane spanned by k and n̂0 and ê2 perpendicular to it. Thus, the wave vector k=k1ê1+k3n̂0 is used to rewrite and summarize Eqs. (A10)(A12) to obtain

Fel=12Vkα=12|n˜α(k)|2(K3k32+Kαk12).
(A13)

It is helpful to remember that n˜α(k) is the Fourier transform of the components of n̂(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

|n˜1(k)|2=VkBTK3k32+K1k12,
(A14)
|n˜2(k)|2=VkBTK3k32+K2k12
(A15)

by analogy with Allen and Frenkel58,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)]

Qαβ=12S(nαnβδαβ)
(A16)

and remember that n̂(r)=(δn1,δn2,1)T [see Eq. (A2)]. Substituting n̂(r) in Eq. (A16) results in Q13(r)=32Sδn1 and Q23(r)=32Sδn2, where the computation of δn1 and δn2 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.

1.
M.
Humar
,
M.
Ravnik
,
S.
Pajk
, and
I.
Muševič
,
Nat. Photonics
3
,
595
(
2009
).
2.
M.
Humar
and
I.
Muševič
,
Opt. Express
18
,
26995
(
2010
).
3.
K.
Peddireddy
,
V. S. R.
Jampani
,
S.
Thutupalli
,
S.
Herminghaus
,
C.
Bahr
, and
I.
Muševič
,
Opt. Express
21
,
30233
(
2013
).
4.
T.
Amann
,
C.
Dold
, and
A.
Kailer
,
Tribol. Int.
65
,
3
(
2013
).
5.
A. M.
Donald
,
A. H.
Windle
, and
S.
Hanna
,
Liquid Crystalline Polymers
(
Cambridge University Press
,
2006
).
6.
A. D.
Rey
,
Soft Matter
6
,
3402
(
2010
).
7.
Y.
Bouligand
, in
Liquid Crystalline Order in Polymers
, edited by
A.
Blumstein
(
Academic Press
,
New York
,
1978
), Chap. 8, pp.
261
297
.
8.
R.
Lipowsky
,
Nature
349
,
475
(
1991
).
9.
A. G.
Petrov
,
Biochim. Biophys. Acta, Biomembr.
1561
,
1
(
2002
).
10.
K.
Luby-Phelps
,
Int. Rev. Cytol.
192
,
189
(
1999
).
11.
J.
Alvarado
,
B. M.
Mulder
, and
G. H.
Koenderink
,
Soft Matter
10
,
2354
(
2014
).
12.
G. T.
Stewart
,
Liq. Cryst.
30
,
541
(
2003
).
13.
S. J.
Woltman
,
G. D.
Jay
, and
G. P.
Crawford
,
Nat. Mater.
6
,
929
(
2007
).
14.
G. T.
Stewart
,
Liq. Cryst.
31
,
443
(
2004
).
15.
A.
Sengupta
,
S.
Herminghaus
, and
C.
Bahr
,
Liq. Cryst. Rev.
2
,
73
(
2014
).
16.
P.-G.
de Gennes
and
J.
Prost
,
The Physics of Liquid Crystals
(
Clarendon
,
Oxford
,
1995
).
17.
M.
Kléman
and
O. D.
Lavrentovich
,
Soft Matter Physics: An Introduction
(
Springer
,
Berlin
,
2007
).
18.
R. S.
Porter
,
E. M.
Barrall
II
, and
J. F.
Johnson
,
J. Chem. Phys.
45
,
1452
(
1966
).
19.
W.
Helfrich
,
Phys. Rev. Lett.
23
,
372
(
1969
).
20.
P. C.
Martin
,
O.
Parodi
, and
P. S.
Pershan
,
Phys. Rev. A
6
,
2401
(
1972
).
21.
P. G.
de Gennes
,
Phys. Fluids
17
,
1645
(
1974
).
22.
L.
Bennett
and
S.
Hess
,
Phys. Rev. E
60
,
5561
(
1999
).
23.
A. J.
Walker
and
I. W.
Stewart
,
Int. J. Eng. Sci.
48
,
1961
(
2010
).
24.
M.
Goulian
and
S. T.
Milner
,
Phys. Rev. Lett.
74
,
1775
(
1995
).
25.
K.
Negita
,
M.
Inoue
, and
S.
Kondo
,
Phys. Rev. E
74
,
051708
(
2006
).
26.
K.
Negita
and
H.
Kaneko
,
Phys. Rev. E
80
,
011705
(
2009
).
27.
S.
Fujii
,
Y.
Ishii
,
S.
Komura
, and
C.-Y. D.
Lu
,
EPL
90
,
64001
(
2010
).
28.
S.
Chatterjee
and
S. L.
Anna
,
Soft Matter
8
,
6698
(
2012
).
29.
C.
Meyer
,
S.
Asnacios
,
C.
Bourgaux
, and
M.
Kleman
,
Rheol. Acta
39
,
223
(
2000
).
30.
D. B.
Liarte
,
M.
Bierbaum
,
M.
Zhang
,
B. D.
Leahy
,
I.
Cohen
, and
J. P.
Sethna
,
Phys. Rev. E
92
,
062511
(
2015
).
31.
A. V.
Kityk
,
M.
Wolff
,
K.
Knorr
,
D.
Morineau
,
R.
Lefort
, and
P.
Huber
,
Phys. Rev. Lett.
101
,
187801
(
2008
).
32.
S.
Calus
,
M.
Busch
,
A. V.
Kityk
,
W.
Piecek
, and
P.
Huber
,
J. Phys. Chem. C
120
,
11727
(
2016
).
33.
J.
Janik
,
R.
Tadmor
, and
J.
Klein
,
Langmuir
13
,
4466
(
1997
).
34.
M.
Ruths
and
S.
Granick
,
Langmuir
16
,
8368
(
2000
).
35.
P.
Huber
,
J. Phys.: Condens. Matter
27
,
103102
(
2015
).
36.
M.
Greschek
,
M.
Melle
, and
M.
Schoen
,
Soft Matter
6
,
1898
(
2010
).
37.
M.
Greschek
and
M.
Schoen
,
Phys. Rev. E
83
,
011704
(
2011
).
38.
S.
Giura
and
M.
Schoen
,
Phys. Rev. E
90
,
022507
(
2014
).
39.
C. G.
Gray
and
K. E.
Gubbins
,
Theory of Molecular Fluids
(
Clarendon
,
Oxford
,
1984
), Vol.
1
.
40.
J. A.
Castellano
,
Mol. Cryst. Liq. Cryst.
94
,
33
(
1983
).
41.
A. A.
Sonin
,
The Surface Physics of Liquid Crystals
(
Gordon and Breach
,
Amsterdam
,
1995
).
42.
K.
Zuhail
,
P.
Sathyanarayana
,
D.
Seč
,
S.
Čopar
,
M.
Škarabot
,
I.
Muševič
, and
S.
Dhara
,
Phys. Rev. E
91
,
030501
(
2015
).
43.
K.
Zuhail
and
S.
Dhara
,
Appl. Phys. Lett.
106
,
211901
(
2015
).
44.
K.
Zuhail
,
S.
Čopar
,
I.
Muševič
, and
S.
Dhara
,
Phys. Rev. E
92
,
052501
(
2015
).
45.
H.
Qi
and
T.
Hegmann
,
J. Mater. Chem.
16
,
4197
(
2006
).
46.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Oxford University Press
,
Oxford
,
1989
).
47.
J. M.
Ilnytskyi
and
M. R.
Wilson
,
Comput. Phys. Commun.
148
,
43
(
2002
).
48.
S.
Nosé
,
J. Chem. Phys.
81
,
511
(
1984
).
49.
S.
Nosé
,
Mol. Phys.
52
,
255
(
1984
).
50.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
51.
R.
Eppenga
and
D.
Frenkel
,
Mol. Phys.
52
,
1303
(
1984
).
52.

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

53.
I. N.
Bronstein
,
K. A.
Semendjajew
,
G.
Musiol
, and
H.
Mühlig
,
Handbuch der Mathematik
(
Harri Deutsch
,
Frankfurt a. M.
,
2000
).
54.
M. G.
Mazza
,
M.
Greschek
,
R.
Valiullin
, and
M.
Schoen
,
Phys. Rev. E
83
,
051704
(
2011
).
55.
S.
Püschel-Schlotthauer
,
T.
Stieger
,
M.
Melle
,
M. G.
Mazza
, and
M.
Schoen
,
Soft Matter
12
,
469
(
2016
).
56.
J.-P.
Hansen
and
I. R.
McDonald
,
Theory of Simple Liquids
, 3rd ed. (
Academic Press
,
London
,
2006
).
57.
M.
Cifelli
and
C.
Veracini
,
Mol. Cryst. Liq. Cryst.
576
,
127
(
2013
).
58.
M. P.
Allen
and
D.
Frenkel
,
Phys. Rev. A
37
,
1813
(
1988
).
59.
M. P.
Allen
and
D.
Frenkel
,
Phys. Rev. A
42
,
3641
(
1990
).
60.
M. P.
Allen
,
M. A.
Warren
,
M. R.
Wilson
,
A.
Sauron
, and
W.
Smith
,
J. Chem. Phys.
105
,
2850
(
1996
).
61.
N.
Madhusudana
and
R.
Pratibha
,
Mol. Cryst. Liq. Cryst.
89
,
249
(
1982
).
62.
P.
De Gennes
,
Solid State Commun.
10
,
753
(
1972
).
63.
L.
Cheung
,
R. B.
Meyer
, and
H.
Gruler
,
Phys. Rev. Lett.
31
,
349
(
1973
).
64.
L.
Leger
,
Phys. Lett. A
44
,
535
(
1973
).
65.
M.
Delaye
,
R.
Ribotta
, and
G.
Durand
,
Phys. Rev. Lett.
31
,
443
(
1973
).
66.
A.
De Vries
,
Mol. Cryst. Liq. Cryst.
11
,
361
(
1970
).
67.
H.-H.
Graf
,
H.
Kneppe
, and
F.
Schneider
,
Mol. Phys.
77
,
521
(
1992
).
68.
F.-J.
Bock
,
H.
Kneppe
, and
F.
Schneider
,
Liq. Cryst.
1
,
239
(
1986
).
69.
B.
Todd
,
D. J.
Evans
, and
P. J.
Daivis
,
Phys. Rev. E
52
,
1627
(
1995
).
70.
K. P.
Travis
,
B.
Todd
, and
D. J.
Evans
,
Phys. Rev. E
55
,
4288
(
1997
).
71.
K. P.
Travis
and
K. E.
Gubbins
,
J. Chem. Phys.
112
,
1984
(
2000
).
72.
T.
Stieger
,
M.
Schoen
, and
M. G.
Mazza
,
J. Chem. Phys.
140
,
054905
(
2014
).
73.
P. M.
Chaikin
and
T. C.
Lubensky
, in
Principles of Condensed Matter Physics
(
Cambridge University Press
,
2000
), Chap. 6.3, pp.
308
316
.
74.
P.
Oswald
and
P.
Pieranski
,
Nematic and Cholesteric Liquid Crystals: Concepts and Physical Properties Illustrated by Experiments
, 1st ed.
Liquid Crystals Book Series
(
CRC Press
,
2005
).