On the elusive saddle–splay and splay–bend elastic constants of nematic liquid crystals

Revignas, Davide; Ferrarini, Alberta

doi:10.1063/5.0153831

The elastic behavior of nematics is commonly described in terms of the three so-called bulk deformation modes, i.e., splay, twist, and bend. However, the elastic free energy contains also other terms, often denoted as saddle–splay and splay–bend, which contribute, for instance, in confined systems. The role of such terms is controversial, partly because of the difficulty of their experimental determination. The saddle–splay (K₂₄) and splay–bend (K₁₃) elastic constants remain elusive also for theories; indeed, even the possibility of obtaining unambiguous microscopic expressions for these quantities has been questioned. Here, within the framework of Onsager theory with Parsons–Lee correction, we obtain microscopic estimates of the deformation free energy density of hard rod nematics in the presence of different director deformations. In the limit of a slowly changing director, these are directly compared with the macroscopic elastic free energy density. Within the same framework, we derive also closed microscopic expressions for all elastic coefficients of rodlike nematics. We find that the saddle–splay constant K₂₄ is larger than both K₁₁ and K₂₂ over a wide range of particle lengths and densities. Moreover, the K₁₃ contribution comes out to be crucial for the consistency of the results obtained from the analysis of the microscopic deformation free energy density calculated for variants of the splay deformation.

I. INTRODUCTION

The Oseen–Frank free energy density of (non-chiral) nematics is expressed in terms of the invariants under nematic symmetry, built from the first-order derivatives of the position (R) dependent director,

\hat{n} (R)

⁠,^1,2

\begin{align} Δ a^{def} (R) & = \frac{1}{2} K_{11} S^{2} + \frac{1}{2} K_{22} T^{2} + \frac{1}{2} K_{33} B^{2} \\ - K_{24} \nabla \cdot (S + B) . \end{align}

(1)

Here,

S = \hat{n} (\nabla \cdot \hat{n})

⁠,

T = \hat{n} \cdot (\nabla \times \hat{n})

⁠, and

B = - (\hat{n} \cdot \nabla) \hat{n} = \hat{n} \times (\nabla \times \hat{n})

represent the splay, twist, and bend deformation modes, respectively, and K₁₁, K₂₂, and K₃₃ are the corresponding elastic constants. The fourth term accounts for the saddle–splay mode, and in the case of a homogeneous degree of order, its integral over a nematic sample can be reduced to a surface integral over the boundary. Hence, it is generally regarded as surface-like elasticity, in contrast to the bulk elasticity described by the other terms. It may be worth mentioning that different definitions of the saddle–splay contribution to the free energy density can be found in the literature. Some authors include a factor of 1/2 in the saddle–splay term in Eq. .^3–5 In other cases, the following alternative form has been employed:^6,7

- \frac{1}{2} (K_{22} + K_{24}^{'}) \nabla \cdot (S + B),

(2)

which compared with Eq. leads to

K_{24}^{'} = 2 K_{24} - K_{22}

⁠. Therefore, attention must be paid when comparing values reported in the literature. Hereafter, data taken from different sources will be converted, when needed, to be consistent with the form of Eq. .

Given its form as a total derivative, the saddle–splay mode does not enter the Euler–Lagrange equations for fixed boundary conditions; however, in principle, it can contribute to the energy through boundary conditions at the surface. It has the effect of renormalizing the surface anchoring, and the difficulty of disentangling the two contributions, together with the negligible role of the saddle–splay term for strong anchoring, are reasons why this mode is often neglected. However, it can have significant implications in complex geometries, such as stripes, cylinders, droplets, or hybrid aligned cells, through the coupling to the local curvature.⁸ The stability of the uniform nematic phase requires that each of the bulk elastic constants is positive, whereas the saddle–splay constant must satisfy the inequalities 0 < K₂₄ < min(K₁₁, K₂₂).⁷ In particular, K₂₂ − K₂₄ < 0 would promote a double twist deformation, as in blue phases.⁹

In a recent reformulation of nematic elasticity,¹⁰ based on a mathematical construction proposed in Ref. 11, the director gradients are decomposed into four modes: in addition to splay, twist, and bend, there is a fourth mode, denoted as Δ, which is directly related to K₂₄. Thus, a form of the free energy density is proposed, with the saddle–splay contribution explicitly regarded as a bulk term.

In addition to the saddle–splay, also another surface-like elastic mode has been considered, which can be expressed as a sum of products of first derivatives of the director and linear terms in its second derivatives.^10,12 It was neglected by Frank,² but it was included by Oseen,¹ Zocher,¹³ and Nehring and Saupe.⁶ This term is generally denoted as (mixed) splay–bend, and its contribution to the free energy density can be expressed as

+ K_{13} \nabla \cdot S,

(3)

where K₁₃ is the corresponding elastic constant. Most of the considerations presented above on the saddle–splay term can be repeated for the splay–bend contribution, with additional questions related to the different geometrical structure of the two modes: the former contains only director derivatives tangential to the surface, whereas in the latter there are also normal derivatives. These are responsible for the unboundedness of the free energy functional, which was a highly debated issue in the last decade of the past century.¹⁴

Bulk elastic constants can be measured by independently exciting the relative modes, through a suitable use of sample geometry, boundary conditions and external fields.¹⁵ On the contrary, surface-like constants are generally evaluated by the analysis of relatively complex director configurations using a certain model. This is one reason why it is difficult to obtain accurate estimates; indeed, very few data are available for thermotropic liquid crystals, and in general, they are affected by a considerable uncertainty. For PCH7 positive, K₁₃, comparable with the bulk constants, was determined in hybrid nematic cells,¹⁶ whereas for 5CB K₁₃ ≈ −0.20K₁₁ was inferred from analysis of the periodicity of the strip domain phase of a hybrid aligned nematic film with azimuthally degenerate boundary conditions.⁴ For the saddle–splay constant K₂₄ some more data, yet rather controversial, are available. K₂₄ of some thermotropic liquid crystals was obtained from the analysis of spontaneously modulated structures in hybrid aligned nematic layers^4,17 or from nuclear magnetic resonance^3,18–20 and optical polarizing microscopy⁵ in systems confined to submicrometer cavities. Despite some discrepancies between the reported data, values comparable in magnitude to the bulk constants were estimated for 5CB, 8CB, and E7. More recently, a new route for measuring K₂₄ was proposed, based on the use of combined topographic and chemical surface patterning to create saddle–splay distortions: K₂₄ ≈ 0.5K₂₂ was estimated for 5CB and for 8CB.²¹ Further support to relatively high values of the saddle–splay constant of 5CB comes from studies of the equilibrium configuration of the nematic director in confined systems: K₂₄ ≈ K₂₂²² was inferred from the formation of doubly twisted configurations in toroidal droplets, whereas K₂₄ ≈ 0.7K₁₁ was estimated from the analysis of the radial to bipolar ordering transition, induced by endotoxin in spherical droplets.²³ Particularly impressive in this respect are the results obtained for chromonics, i.e., water solutions of plank-like aromatic molecules with charged substituents in the periphery, which self-assemble into columnar stacks and beyond a certain concentration form liquid crystal phases.²⁴ Disodium cromoglycate (DSCG) and Sunset Yellow (SSY), two well-studied representatives of this class of liquid crystals, were found to adopt an escaped twist configuration, typical of blue phases,⁹ when confined in cylindrical capillaries with degenerate planar boundary conditions^25–27 and in cylindrical shells.²⁸ This was surprising because the molecules are achiral and could only be explained by invoking a stabilization due to the saddle–splay contribution to the free energy. Particularly high K₂₄ values were inferred from the analysis of the director configurations in these systems: K₂₄/K₂₂ = 27.5²⁵ and 3.25^29,30 for SSY, K₂₄/K₂₂ = 7.5^27,28 for DSCG.

Molecular theories and simulations, which in principle could provide independent insights on the relation between the various elastic constants and their dependence on the molecular structure and interactions, have devoted scarce attention to surface-like elasticity. Within molecular theories, expressions for elastic constants of nematics are obtained from the free energy in the presence of director deformations, which are generally analyzed through Taylor expansion of the director as a function of position. Early estimates of both the bulk and the surface-like elastic constants were provided by Nehring and Saupe, based on anisotropic dispersion interactions (with implicitly isotropic short-range repulsions),³¹ and the relationships K₁₁:K₂₂:K₃₃:K₁₃ = 5:11:5:−6 were obtained, together with very small positive K₂₄. Remarkably, the very possibility of unambiguously defining K₁₃ and K₂₄ from a microscopic theory was questioned.³² This objection was discussed in a subsequent study,³³ where the surface-like elastic constants of a Gay–Berne fluid were evaluated using density functional theory. K₁₃ and K₂₄ close to K₂₂/2 were predicted, and the discrepancy from previous results^31,34 was ascribed to the fact that in previous studies the pair potential favored end-to-end over side-by-side alignment. Moreover, the results could be affected by the level of approximation adopted in the theory, first of all the truncation at the second rank orientational order parameter in the series expansion of the single particle density, which leads to K₃₃ = K₁₁. In Ref. 35, dealing again with a Gay–Berne system, order parameters up to fourth rank were retained in the spherical harmonics expansion of the direct correlation function, which was calculated from molecular dynamics data for the pair correlation function of the nematic phase by solving the Ornstein–Zernicke equation, rather than from a simple isotropic liquid approximation. At variance with the previous study, K₂₄ close to K₂₂ and small negative K₁₃ were obtained, the latter almost an order of magnitude smaller, in absolute value, than the bulk constants. The surface-like elastic constants were found to be more sensitive to details of the pair correlation function than the bulk ones. Finally, in Ref. 36, an atomistic model that can be seen as an extension of the Maier–Saupe potential was used for some typical thermotropic systems. The prediction was K₂₄ > K₂₂ for PAA, and K₂₄ ≈ K₂₂ for 5CB and 8CB, whereas K₁₃ was found to be small and positive for PAA, small and negative for 5CB and even more negative for 8CB. The strong dependence on the molecular structure suggests possible roles of the molecular stiffness and biaxiality, which increase on going from 8CB to PAA. Recently, numerical estimates of the saddle–splay elastic constant were provided by molecular dynamics simulations. For Gay–Berne particles, K₂₄ ≈ K₂₂ was obtained,³⁷ whereas for an atomistic model of 5CB, K₂₄ between 1.5 and 2.5 pN, smaller than K₂₂, was found.³⁸

The scenario is confusing: the magnitude of the surface-like elastic constants remains unclear and there is no insight into their dependence on the thermodynamic parameters and the microscopic structure of materials. This is quite different from the case of the bulk elastic constants, whose behavior is rather well understood, not only in conventional nematics, which typically exhibit K₂₂ < K₁₁ < K₃₃, but also in less conventional systems.^39–41 Recently,⁴² we addressed the question of unambiguous determination of the saddle–splay elastic constant in the framework of Onsager–Straley theory^43,44 for hard rodlike particles. K₂₄ and the three bulk elastic constants were calculated from the profile of the deformation free energy as a function of the deformation wavenumber, for a double twist distortion and for the standard splay, twist, and bend modes, respectively: K₂₄ > K₂₂ was found, which implies an intrinsic tendency to twist.

Here, we extend the approach proposed in Ref. 42 to other deformation modes, which allow us to obtain estimates also for the splay–bend constant K₁₃. Moreover, we derive explicit expressions for all the elastic constants using a second-order Taylor expansion of the singlet orientational distribution function (ODF). We then report numerical estimates of the full set of elastic constants for tangent hard sphere rods (THSRs) with different aspect ratios and over a wide range of densities.

This paper is organized as follows. Section II is devoted to the theoretical development: first the deformation modes are defined using continuum theory; then, a microscopic expression for the deformation free energy density is obtained, and finally, expressions for elastic constants are derived. Section III is dedicated to the numerical implementation of the theoretical expressions. In Sec. IV, we report the results of the different calculations and finally, in Sec. V, we summarize the conclusions of our study.

II. THEORY

A. Macroscopic deformation analysis

We start our analysis by introducing a set of vector functions of the position R (expressed with respect to a fixed laboratory frame OXYZ) that describes different deformations of a uniform director field. The explicit forms of these deformations are listed in Table I. All deformations share three features: (i) their definition involves a parameter, indicated as q, which is a measure of the degree of deformation; (ii) for q = 0, a uniform director field

{\hat{n}}_{u} ‖ \hat{Z}

is recovered; and (iii) the director at the origin is in any case parallel to the laboratory Z axis,

{\hat{n}}_{d} (0, q) = \hat{Z}

⁠. Using these expressions, the elastic free energy density, sum of Eqs. and , at the origin of the reference frame, R = 0, has a parabolic dependence on q and the parabolic coefficients are given by linear combinations,

C_{d}

⁠, of elastic constants

Δ a_{d}^{def} (0, q) = q^{2} C_{d} (K_{11}, K_{22}, K_{33}, K_{24}, K_{13}),

(4)

where the subscript “d” stands for one of the director fields listed in Table I; the corresponding explicit expressions of

C_{d}

are also listed there. Three functions (splay, twist, and bend) are linked each to a single bulk elastic constant. Three additional functions were proposed in Ref. 10; one of them, Δ [see Fig. 1(a)], involves only K₂₄, whereas the other two, double twist and double splay [see Fig. 1(b)], are related to linear combinations of K₂₄, K₁₁, and K₂₂. Here, we define a generalized version of double splay in which the additional parameter w has been introduced in order to control the degree of involvement of the K₁₃ constant.

TABLE I.

Deformed director field functions ${\hat{n}}_{d} (R, q)$ employed in this work and corresponding combinations of elastic constants $(C_{d})$ ⁠.

Deformation name (d)	${\hat{n}}_{d} (R, q)$	$C_{d} = Δ a_{d}^{def} (0, q) / q^{2}$
Splay (S)	$[q R_{X}, 0, 1 + q R_{Z}] / \sqrt{q^{2} R_{X}^{2} + {(1 + q R_{Z})}^{2}}$	K₁₁/2
Twist (T)	[sin(qR_Y), 0, cos(qR_Y)]	K₂₂/2
Bend (B)	$[q R_{Z}, 0, 1 - q R_{X}] / \sqrt{q^{2} R_{Z}^{2} + {(1 - q R_{X})}^{2}}$	K₃₃/2
Δ deformation (Δ)	$[q R_{X}, - q R_{Y}, 1] / \sqrt{q^{2} (R_{X}^{2} + R_{Y}^{2}) + 1}$	2K₂₄
Double twist (DT)	$[- q R_{Y}, q R_{X}, 1] / \sqrt{q^{2} (R_{X}^{2} + R_{Y}^{2}) + 1}$	2(K₂₂ − K₂₄)
Generalized	$[q R_{X}, q R_{Y}, 1 + q w R_{Z}] / \sqrt{q^{2} (R_{X}^{2} + R_{Y}^{2}) + {(1 + q w R_{Z})}^{2}}$	2[K₁₁ + (2 − w)K₁₃ − K₂₄]a
double splay (DS)

Deformation name (d)	${\hat{n}}_{d} (R, q)$	$C_{d} = Δ a_{d}^{def} (0, q) / q^{2}$
Splay (S)	$[q R_{X}, 0, 1 + q R_{Z}] / \sqrt{q^{2} R_{X}^{2} + {(1 + q R_{Z})}^{2}}$	K₁₁/2
Twist (T)	[sin(qR_Y), 0, cos(qR_Y)]	K₂₂/2
Bend (B)	$[q R_{Z}, 0, 1 - q R_{X}] / \sqrt{q^{2} R_{Z}^{2} + {(1 - q R_{X})}^{2}}$	K₃₃/2
Δ deformation (Δ)	$[q R_{X}, - q R_{Y}, 1] / \sqrt{q^{2} (R_{X}^{2} + R_{Y}^{2}) + 1}$	2K₂₄
Double twist (DT)	$[- q R_{Y}, q R_{X}, 1] / \sqrt{q^{2} (R_{X}^{2} + R_{Y}^{2}) + 1}$	2(K₂₂ − K₂₄)
Generalized	$[q R_{X}, q R_{Y}, 1 + q w R_{Z}] / \sqrt{q^{2} (R_{X}^{2} + R_{Y}^{2}) + {(1 + q w R_{Z})}^{2}}$	2[K₁₁ + (2 − w)K₁₃ − K₂₄]a
double splay (DS)

^a

The combination $C_{DS}$ for the case w = 1 reported here differs from the one in Ref. 10 because there the K₁₃ term was not included in the macroscopic free energy density expression.

FIG. 1.

View large Download slide

Schematic illustration of the director deformations (a) Δ and (b) generalized double splay, defined in Table I, in the plane Z = 0. In this plane, the splayed director field is independent of the value of the w parameter.

B. Microscopic free energy density of deformation

In this section, we obtain a microscopic, particle-based equivalent of Eq. and we describe a method to estimate the elastic constants. In the framework of Onsager theory,⁴³ the Helmholtz free energy of a system of uniaxial particles is expressed as functional of the single particle ODF, i.e., the probability density of finding a particle with its center of mass (c.m.) at position R and with orientation

\hat{u}

⁠. For the general case of a nonuniform director field

{\hat{n}}_{d}

⁠, we will use the following ansatz for the form of the ODF:

f_{d} (\hat{u}, R, q) = \frac{\exp {c {[\hat{u} \cdot {\hat{n}}_{d} (R, q)]}^{2}}}{Q},

(5)

where c is a parameter that quantifies the orienting strength and

Q = \int d \hat{u} \exp {c {[\hat{u} \cdot {\hat{n}}_{d} (R, q)]}^{2}}

is the orientational partition function. Equation is the simplest functional form compatible with the uniaxial symmetry of the particle and the local nematic symmetry. Inclusion of higher rank terms may affect the numerical results at high ordering, but would be unessential for the scope of this work. The dependence on the scalar product

[\hat{u} \cdot {\hat{n}}_{d}]

implies that, at fixed c, the ODF does not change its profile when a deformation is present, but it simply rotates to match the local director.

In the case of a uniform director field

{\hat{n}}_{u} ‖ \hat{Z}

⁠, i.e., for an undeformed nematic system, the ODF takes the simple form

f_{u} (\hat{u}) = \frac{\exp {c u_{Z}^{2}}}{Q} .

(6)

The Helmholtz free energy of the system is given by the sum of the ideal and the excess part. The former accounts for the decrease in rotational entropy due to nematic ordering and reduces to the isotropic ideal gas term when the ODF is constant. Its expression per unit volume for an undeformed system is

\frac{a_{u}^{id}}{k_{B} T} = \frac{a_{u}^{iso}}{k_{B} T} + ρ \int d \hat{u} f_{u} (\hat{u}) \ln [4 π f_{u} (\hat{u})],

(7)

where a^iso is the isotropic ideal gas free energy density, k_B is Boltzmann’s constant, T is the temperature of the system, and ρ = N/V is the density of particles, N being their total number enclosed in a volume V. The ideal free energy density in the presence of a deformation,

a_{d}^{id}

⁠, is obtained by substituting the general form of the ODF [Eq. ] in Eq. ; however, since the latter involves the integration over all possible orientations and the ODF is evaluated at a single point in space,

a_{d}^{id} = a_{u}^{id}

is obtained.

The excess contribution to the Helmholtz free energy, which accounts for the interactions between particles, is expressed according to second virial Onsager theory with modified⁴⁵ Parsons–Lee^46,47 correction. In particular, the excess free energy density in R = 0 for the deformations defined in Table I reads

\begin{align} \frac{a_{d}^{ex} (0, q)}{k_{B} T} & = - \frac{η ρ^{2}}{2} \int d {\hat{u}}_{A} \int d {\hat{u}}_{B} \int d R_{A B} \\ \times f_{d} ({\hat{u}}_{A}, 0, q) f_{d} ({\hat{u}}_{B}, R_{A B}, q) e_{A B} ({\hat{u}}_{A}, {\hat{u}}_{B}, R_{A B}) \end{align}

(8)

where A and B are labels for two particles at positions R_A = 0 and R_B, respectively, and R_AB = R_B − R_A is the vector joining their c.m.s; e_AB is the Mayer function,⁴⁸ which in turn is defined as e_AB = exp{−U_AB/k_BT} − 1 with U_AB being the interaction potential energy between the two particles. The triple integral represents the volume excluded to particle B by particle A, averaged over all their possible orientations. Finally, η is the modified Parsons–Lee^46,47 factor

η = \frac{1 - \frac{3}{4} ρ v_{eff}}{{(1 - ρ v_{eff})}^{2}},

(9)

which corrects for neglecting higher order virial terms; v_eff is an effective volume, defined as the portion of space in which any part of a particle cannot enter due to the presence of another particle.⁴⁵ In the undeformed state, the ODFs f_d in Eq. are replaced by the corresponding position independent ODFs, f_u, defined in Eq. .

As explained in Appendix A, the difference between the excess free energy density in the presence and in the absence of director distortions provides a microscopic expression for the (macroscopic) deformation free energy density [Eqs. and ]

Δ a_{d}^{def} (0, q) = a_{d}^{ex} (0, q) - a_{u}^{ex} .

(10)

Therefore, for a given form of the director field function

{\hat{n}}_{d}

⁠, Eq. can be used to calculate the deformation free energy density as a function of the parameter q, with the only assumption of constant orienting strength c. This in turn corresponds to a fixed value of the nematic order parameter

S = (3 {⟨ u_{Z}^{2} ⟩}_{u} - 1) / 2

⁠, where ⟨⋯⟩_u indicates the orientational average over the undeformed ODF defined in Eq. . From the profile of

Δ a_{d}^{def} (0, q)

⁠, one can obtain insight on the adequacy of the first or second-order form of the elastic free energy density and the possible role of higher order contributions,⁴⁹ whereas from the curvature of

Δ a_{d}^{def} (0, q)

for q → 0, we can determine the elastic constants (see the expressions of

C_{d}

in Table I).

C. Explicit expressions for the elastic constants

A second route for the determination of the elastic constants, closer to the conventional one, is based on the Taylor expansion around q = 0 of the ODFs appearing in Eq. . If we introduce the following shorthand notations:

\begin{gathered} f_{d}^{A} = f_{d} ({\hat{u}}_{A}, 0, q) = f_{u} ({\hat{u}}_{A}) = f_{u}^{A} \\ f_{d}^{B} = f_{d} ({\hat{u}}_{B}, R_{A B}, q), f_{u}^{B} = f_{u} ({\hat{u}}_{B}) \\ e_{A B} = e_{A B} ({\hat{u}}_{A}, {\hat{u}}_{B}, R_{A B}), \end{gathered}

we can approximate the product of ODFs in the integral of Eq. as

\begin{align} f_{d}^{A} f_{d}^{B} & ≃ f_{u}^{A} f_{u}^{B} + q \{[\partial_{q} f_{d}^{A} |_{q = 0}] f_{u}^{B} \\ + [\partial_{q} f_{d}^{B} |_{q = 0}] f_{u}^{A}\} + \frac{q^{2}}{2} \{[\partial_{q} \partial_{q} f_{d}^{A} |_{q = 0}] f_{u}^{B} \\ + [\partial_{q} \partial_{q} f_{d}^{B} |_{q = 0}] f_{u}^{A} + 2 [\partial_{q} f_{d}^{A} |_{q = 0}] [\partial_{q} f_{d}^{B} |_{q = 0}]\}, \end{align}

(11)

where the symbol ∂_q indicates differentiation with respect to q. If we introduce Eq. into the expression of the microscopic free energy density [Eq. ], we obtain the following correspondence between the combinations of elastic constants

C_{d}

and the second-order term of the expansion:

\begin{align} \frac{C_{d}}{k_{B} T} & = - \frac{η ρ^{2}}{4} \int d {\hat{u}}_{A} \int d {\hat{u}}_{B} \int d R_{A B} e_{A B} \\ \times \{[\partial_{q} \partial_{q} f_{d}^{A} |_{q = 0}] f_{u}^{B} + [\partial_{q} \partial_{q} f_{d}^{B} |_{q = 0}] f_{u}^{A} \\ + 2 [\partial_{q} f_{d}^{A} |_{q = 0}] [\partial_{q} f_{d}^{B} |_{q = 0}]\} . \end{align}

(12)

In order to obtain expressions for the elastic constants, we must make the ODF derivatives appearing in Eq. explicit. The first derivative of the ODF evaluated at q = 0 has the general form

[\partial_{q} f_{d} |_{q = 0}] = f_{u} [\partial_{q} g_{d} |_{q = 0}],

(13)

where we introduced the shorthand notation

g_{d} = g_{d} (\hat{u}, R, q) = c {[\hat{u} \cdot {\hat{n}}_{d} (R, q)]}^{2} .

(14)

Similarly, the second derivative of the ODF evaluated at q = 0 has the general form

[\partial_{q} \partial_{q} f_{d} |_{q = 0}] = f_{u} \{{[\partial_{q} g_{d} |_{q = 0}]}^{2} + [\partial_{q} \partial_{q} g_{d} |_{q = 0}]\} .

(15)

The explicit expressions of the first and the second derivatives of

g_{d} (\hat{u}, R, q)

evaluated at q = 0 are collected in Table II. We note that, since all terms involve at least one Cartesian component of R and particle A is placed at the origin, R_A = 0, only the second term between braces in Eq. survives, so we can write

\frac{C_{d}}{k_{B} T} = - \frac{η ρ^{2}}{4} ⟨\int d R_{A B} e_{A B} \{{[\partial_{q} g_{d}^{B} |_{q = 0}]}^{2} + [\partial_{q} \partial_{q} g_{d}^{B} |_{q = 0}]\}⟩,

(16)

where we have introduced the angular brackets to denote a double orientational average

⟨\dots⟩ = \int d {\hat{u}}_{A} \int d {\hat{u}}_{B} f_{u}^{A} f_{u}^{B} [\dots] .

(17)

For example, if we use the splay deformation

{\hat{n}}_{S}

from Table I, we obtain an expression for

K_{11} = 2 C_{S}

⁠,

\begin{align} \frac{K_{11}}{k_{B} T} & = - \frac{η ρ^{2}}{2} ⟨\int d R e_{A B} \{4 c^{2} R_{X}^{2} u_{B X}^{2} u_{B Z}^{2} \\ + 2 c R_{X}^{2} (u_{B X}^{2} - u_{B Z}^{2}) - 4 c R_{X} R_{Z} u_{B X} u_{B Z}\}⟩, \end{align}

(18)

where we dropped the “AB” subscript from R_AB for simplicity. If we now introduce the following notation for the opposite of the spatial integral of the Mayer function multiplied by a generic factor:

E_{[\dots]} ({\hat{u}}_{A}, {\hat{u}}_{B}) = - \int d R e_{A B} ({\hat{u}}_{A}, {\hat{u}}_{B}, R) [\dots],

(19)

we can further simplify Eq. obtaining

\begin{aligned} \frac{K_{11}}{k_{B} T} & = η ρ^{2} c [2 c ⟨ E_{R_{X}^{2}} u_{B X}^{2} u_{B Z}^{2} ⟩ + ⟨ E_{R_{X}^{2}} (u_{B X}^{2} - u_{B Z}^{2}) ⟩ \\ - 2 ⟨ E_{R_{X} R_{Z}} u_{B X} u_{B Z} ⟩] . \end{aligned}

(20)

Similarly, from the twist and the bend deformations, we obtain the expressions for K₂₂ and K₃₃,

\frac{K_{22}}{k_{B} T} = η ρ^{2} c [2 c 〈 E_{R_{Y}^{2}} u_{B X}^{2} u_{B Z}^{2} 〉 + 〈 E_{R_{Y}^{2}} (u_{B X}^{2} - u_{B Z}^{2}) 〉],

(21)

\begin{align} \frac{K_{33}}{k_{B} T} & = η ρ^{2} c [2 c ⟨ E_{R_{Z}^{2}} u_{B X}^{2} u_{B Z}^{2} ⟩ + ⟨ E_{R_{Z}^{2}} (u_{B X}^{2} - u_{B Z}^{2}) ⟩ \\ + 2 ⟨ E_{R_{X} R_{Z}} u_{B X} u_{B Z} ⟩] . \end{align}

(22)

Using the Δ deformation, we arrive at the following expression for K₂₄:

\begin{align} \frac{K_{24}}{k_{B} T} & = \frac{η ρ^{2} c}{4} \{2 c [⟨E_{R_{X}^{2}} u_{B X}^{2} u_{B Z}^{2}⟩ + ⟨E_{R_{Y}^{2}} u_{B Y}^{2} u_{B Z}^{2}⟩ + \\ - 2 ⟨E_{R_{X} R_{Y}} u_{B X} u_{B Y} u_{B Z}^{2}⟩] + ⟨E_{R_{X}^{2}} (u_{B X}^{2} - u_{B Z}^{2})⟩ \\ + ⟨E_{R_{Y}^{2}} (u_{B Y}^{2} - u_{B Z}^{2})⟩ - 2 ⟨E_{R_{X} R_{Y}} u_{B X} u_{B Y}⟩\} . \end{align}

(23)

By exploiting Eqs. and (see Appendix B), the last expression can be further simplified to

\begin{align} \frac{K_{24}}{k_{B} T} & = \frac{η ρ^{2} c}{2} [2 c ⟨E_{R_{X}^{2}} u_{B X}^{2} u_{B Z}^{2}⟩ - 2 c ⟨E_{R_{X} R_{Y}} u_{B X} u_{B Y} u_{B Z}^{2}⟩ \\ + ⟨E_{R_{X}^{2}} (u_{B X}^{2} - u_{B Z}^{2})⟩ - ⟨E_{R_{X} R_{Y}} u_{B X} u_{B Y}⟩] . \end{align}

(24)

From the w-dependent quantities in the generalized double splay deformation, we obtain the following expression for K₁₃:

\frac{K_{13}}{k_{B} T} = \frac{η ρ^{2} c}{2} [⟨E_{R_{X} R_{Z}} u_{B X} u_{B Z}⟩ + ⟨E_{R_{Y} R_{Z}} u_{B Y} u_{B Z}⟩],

(25)

which, by making use of Eq. , can be further simplified to

\frac{K_{13}}{k_{B} T} = η ρ^{2} c ⟨E_{R_{X} R_{Z}} u_{B X} u_{B Z}⟩ .

(26)

Finally, using Eqs. –, the Nehring–Saupe⁶ relationship between K₁₃ and other elastic constants is recovered,

K_{13} = 2 K_{24} - \frac{1}{2} (K_{11} + K_{22}) .

(27)

TABLE II.

First and second derivatives with respect to q of the argument g_d of the exponential function of the ODF for the deformations defined in Table I.

d	$\frac{\partial_{q} {[g_{d} (\hat{u}, R, q)]}_{q = 0}}{2 c}$	$\frac{\partial_{q} \partial_{q} {[g_{d} (\hat{u}, R, q)]}_{q = 0}}{2 c}$
S	R_Xu_Xu_Z	$R_{X}^{2} (u_{X}^{2} - u_{Z}^{2}) - 2 R_{X} R_{Z} u_{X} u_{Z}$
T	R_Yu_Xu_Z	$R_{Y}^{2} (u_{X}^{2} - u_{Z}^{2})$
B	R_Zu_Xu_Z	$R_{Z}^{2} (u_{X}^{2} - u_{Z}^{2}) + 2 R_{X} R_{Z} u_{X} u_{Z}$
Δ	R_Xu_Xu_Z − R_Yu_Yu_Z	$R_{X}^{2} (u_{X}^{2} - u_{Z}^{2}) + R_{Y}^{2} (u_{Y}^{2} - u_{Z}^{2}) - 2 R_{X} R_{Y} u_{X} u_{Y}$
DT	R_Xu_Yu_Z − R_Yu_Xu_Z	$R_{X}^{2} (u_{Y}^{2} - u_{Z}^{2}) + R_{Y}^{2} (u_{X}^{2} - u_{Z}^{2}) - 2 R_{X} R_{Y} u_{X} u_{Y}$
DS	R_Xu_Xu_Z + R_Yu_Yu_Z	$R_{X}^{2} (u_{X}^{2} - u_{Z}^{2}) + R_{Y}^{2} (u_{Y}^{2} - u_{Z}^{2}) + 2 R_{X} R_{Y} u_{X} u_{Y}$
		$- 2 w (R_{X} R_{Z} u_{X} u_{Z} + R_{Y} R_{Z} u_{Y} u_{Z})$

III. NUMERICAL IMPLEMENTATION

A. Model details

In this work, we focused on the simple particle model of THSRs: the particles are straight chains of M = 6, 12, and 24 spheres of diameter σ. The potential energy between two rods A and B is defined as

U_{A B} ({\hat{u}}_{A}, {\hat{u}}_{B}, R_{A B}) = \{\begin{cases} \infty & \begin{matrix} if at least one sphere of A overlaps \\ with one or more spheres of B \end{matrix} \\ 0 & otherwise, \end{cases}

(28)

where

{\hat{u}}_{A}

and

{\hat{u}}_{B}

are unit vectors parallel to the C_∞ axes of the rods and R_AB is the vector joining their centers of mass. The system is athermal, i.e., T has a trivial effect and the thermodynamic state is determined only by the number density ρ.

The geometrical volume v₀ of a particle, used to calculate the packing fraction ϕ = ρv₀ of the system, is simply the sum of the volumes of the hard spheres: v₀ = Mπσ³/6. The effective volume, employed in the evaluation of the Parsons–Lee coefficient η [Eq. (9)] was calculated as the volume enclosed by the tessellated surface defined by a sphere of radius equal to 0.5 σ rolling over the particle;⁵⁰ a density of vertices on the surface equal to 1000 σ⁻² was adopted. The results for the three particle types are reported in Table III.

TABLE III.

Geometric and thermodynamic properties of the THSRs studied in this work: number of hard spheres M, geometric (v₀), and effective (v_eff) volume, packing fractions of the isotropic, ϕ_I, and nematic phase, ϕ_N, at the coexistence.

M	v₀/σ³	v_eff/σ³	ϕ_I	ϕ_N
6	3.14	3.50	0.337	0.354
12	6.28	7.08	0.186	0.207
24	12.57	14.24	0.096	0.113

B. Calculation of integrals

In this work, we are concerned with the calculation of two main type of integrals: free energy densities in presence of a certain deformation [Eq. (8)] and generalized excluded volumes in the form of Eq. (19) orientationally averaged over undeformed ODFs, as required, for instance, in Eq. (20).

For the first type of integral, we adopt a simple Monte Carlo scheme consisting of (i) generation of uniformly distributed random orientations for particles A and B over the unitary sphere and of a uniformly distributed random position for the center of mass of particle B inside a sphere of radius equal to Mσ (particle A’s center of mass is fixed at the origin); (ii) check for overlap between the two particles, i.e., check if the Mayer function e_AB is different from zero; (iii) if the particles do overlap, evaluate the two ODFs inside the integral. The final values of $a_{d}^{ex} (0, q)$ are obtained as simple averages over all the evaluations performed in the previous steps. For all the calculations reported in Figs. 3 and 4, the number of steps employed to obtain convergence is N_MC = 10¹².

Generalized excluded volume integrals over R in the form of Eq. (19) are calculated in spherical coordinates, R = [R, Θ, Φ]. Integration is performed analytically on the radius R and by numerical quadrature over the polar angle Θ and the azimuthal angle Φ. This is repeated on a 4-indexes discrete grid of angles {α_A, β_A, α_B, β_B} defining orientations of the two uniaxial particles, ${\hat{u}}_{A}$ and ${\hat{u}}_{B}$ ⁠. The results are then stored into matrices that can be efficiently employed for the evaluation of orientational averages in the form of Eq. (17) at different values of the orienting strength parameter c. For all polar angles (Θ, β_A, β_B), the Gauss–Legendre quadrature scheme is employed, whereas a Gauss–Chebyshev quadrature is used for azimuthal angles (Φ, α_A, α_B).⁵¹

C. Isotropic–nematic coexistence and equilibrium order parameter

At a certain thermodynamic condition, i.e., for a fixed value of the number density ρ, the equilibrium ODF is obtained by minimizing the total free energy density of the undeformed system, $a_{u} = a_{u}^{ex} + a_{u}^{id}$ ⁠. According to our ansatz, the ODF [Eq. (6)] depends on a single parameter, c, which has to be varied in order to obtain the minimum a_u for a certain ρ value. In this way, a unique relation between the orienting strength c and the number density ρ is established, c = c(ρ). Practically, this is obtained by applying the Nelder–Mead simplex algorithm⁵² implemented in Matlab⁵³ on the sum of Eqs. (7) and (8) for the undeformed system, whose integrals are evaluated with the quadrature procedure outlined in Subsection III B.

Using the relation c(ρ), it is possible to calculate the pressure P and the chemical potential μ of the nematic phase by differentiating the free energy density with respect to the volume V and the number of particles N of the system, respectively,

\begin{matrix} \frac{P (ρ)}{k_{B} T} & = & - \frac{1}{k_{B} T} {(\frac{\partial (V a_{u})}{\partial V})}_{N, T} \\ = & ρ + \frac{ρ^{2}}{2} ⟨ E_{1} ⟩ [\frac{ρ v_{eff} (5 - 3 ρ v_{eff})}{4 {(1 - ρ v_{eff})}^{3}} + η], \end{matrix}

(29)

\begin{matrix} \frac{μ (ρ)}{k_{B} T} & = & \frac{1}{k_{B} T} {(\frac{\partial (V a_{u})}{\partial N})}_{V, T} \\ = & \ln (ρ σ^{3}) + \int d \hat{u} f_{u} (\hat{u}) \ln f_{u} (\hat{u}) \\ + \frac{ρ}{2} ⟨ E_{1} ⟩ [\frac{ρ v_{eff} (5 - 3 ρ v_{eff})}{4 {(1 - ρ v_{eff})}^{3}} + 2 η] + const, \end{matrix}

(30)

where

E_{1} ({\hat{u}}_{A}, {\hat{u}}_{B})

⁠, calculated according to Eq. , is simply the excluded volume calculated for a given orientation of the two particles. The densities of the isotropic ρ_iso and the nematic ρ_nem phases in the coexistence regime are obtained from the solution of the following system of equations:

\{\begin{cases} P (ρ_{nem}) = P_{iso} (ρ_{iso}), \\ μ (ρ_{nem}) = μ_{iso} (ρ_{iso}), \end{cases}

(31)

where P_iso(ρ) and μ_iso(ρ) are the pressure and the chemical potential of the coexisting isotropic phase, obtained by enforcing c = 0 in the ODF. The packing fractions of the coexisting nematic and isotropic phases calculated for the three investigated systems are reported in Table III.

D. Calculations for the contour plots of the average excluded volume

The spatially resolved contributions to the average excluded volume, ⟨E₁⟩, shown in the plots of Fig. 5 for deformations of cylindrical symmetry are calculated as follows: (i) the XZ plane is discretized into a square grid of mesh size b; (ii) the orientational average of the Mayer function e_AB is evaluated in each cell by sampling the orientations of the particles A (fixed at the origin, see Fig. 2) and B (with its c.m. at the center of the current cell) from their ODFs [Eq. (5)]; (iii) each cell value is scaled by a factor equal to the distance from the Z axis so that the total excluded volume would be obtained by simply summing over all cells and multiplying by b²π. Calculations are performed for fixed values of the c and q parameters (q = 0 for the undeformed system).

FIG. 2.

View large Download slide

Sketch of the pair configurations considered to calculate contour plots of the average excluded volume: particle A (in yellow) has its c.m. fixed at the origin whereas particle B (in red) moves in the XZ plane. When the center of mass of particle B is outside the circle of radius Mσ, there cannot be any overlap between the two particles. The circle encloses the positions of the c.m. of B in which overlapping with particle A may occur.

IV. RESULTS

A. Microscopic free energy density for different deformations

Figures 3 and 4(a) show the deformation free energy density [Eq. (10)] as a function of the wavenumber q, calculated for different distortions of the nematic phase of THSRs with M = 12, at volume fraction ϕ = 0.282. In all cases the calculated points were fitted (R² > 0.99) to a parabolic function and from the coefficients the bulk and surface-like elastic constants were evaluated (see Table I). Figure 3 shows $Δ a_{d}^{def} (0, q)$ for splay, twist, bend, Δ, and double twist deformations, whose curvatures provide K₁₁ = 1.54 k_BTσ⁻¹, K₂₂ = 0.60 k_BTσ⁻¹, K₃₃ = 20.18 k_BTσ⁻¹, and K₂₄ = 2.17 k_BTσ⁻¹.

FIG. 3.

View large Download slide

Deformation free energy density [Eq. (10)] in the presence of splay, twist, bend, Δ, and double twist distortions, as a function of the deformation wavenumber q, calculated for THSRs with M = 12 at number density ϕ = 0.282 (symbols). The solid lines show the results from fitting to Eq. (4); the corresponding parabolic coefficients $C_{d}$ are reported in the legend.

FIG. 4.

View large Download slide

(a) Deformation free energy density [Eq. (10)] as a function of the wavenumber q, for generalized double splay director fields (see Table I), corresponding to values of the w parameter ranging from −1 to 2.5 in steps of 0.5 (color from brown to green). The results, obtained for THSRs with M = 12 at number density ϕ = 0.282, are shown as star symbols, while the lines are the results from fitting to Eq. (4). The dotted and the dashed lines correspond to the cases w = 1 (double splay defined in Ref. 10) and w = 2, respectively. (b) Estimated parabolic coefficients $C_{DS}$ from previous fitting (black crosses) as a function of the w parameter. The solid red line in the same plot corresponds to their linear fit according to the expression in Table I; from its slope an estimated value of K₁₃ ≃ 3.29 k_BTσ⁻¹ is obtained.

Interestingly, K₂₄ is larger than K₂₂, which is confirmed by direct calculation of the difference K₂₂ − K₂₄ from the curvature of the deformation free energy density for double twist deformation.

Figure 4(a) shows the deformation free energy as a function of the wavenumber q for generalized splayed director fields, corresponding to values of w (see Table I) ranging from −1 to 2.5 in steps of 0.5. We can see that the deformations have different costs and, in particular, we can observe a change in curvature between w = 1.5 and w = 2. The coefficients $C_{D S} (w)$ obtained from the parabolic fitting in q of these data are reported in Fig. 4(b). From the linear fitting of such coefficients as a function of w according to the expression reported in Table I, we obtain K₁₃ = 3.29 k_BTσ⁻¹. Thus, the saddle–splay constant turns out to be higher than all the other elastic coefficients but K₃₃. Even more important is that the inclusion of the splay–bend term in the elastic free energy, which is ignored by first-order theories, is needed to consistently distinguish the results of microscopic calculations for different variants of the generalized double splay.

A final comment is required by the negative curvatures observed in Fig. 3 for double twist and in Fig. 4(a) for generalized double splay with w = 2 (dashed line), which derive from K₂₄ > K₂₂ and K₂₄ > K₁₁, respectively (in the latter case there is no K₁₃ contribution). This points to a local instability of uniform uniaxial order w.r.t such deformations. Given the form of the interparticle potential [Eq. (28)], the phenomenon has an entropic origin: it is driven by the decrease of excluded volume between a pair of rods fluctuating around a local director, in the presence of these deformations. This is illustrated by the contour plots in Fig. 5, which show the point-by-point difference between the average excluded volume for two THSRs in the presence of a double twist (a) or a generalized double splay with w = 2 (b) director field w.r.t. the same quantity but in a uniform director. In the case of double twist, there is a decrease of average excluded volume for any relative position of the particles, with a more pronounced effect when their c.m.s are close to each other. On the contrary, for generalized double splay, the overall effect comes from the competition between two counteracting tendencies: when the moving particle is in the region where the director field lines diverge, the probability of overlap is higher than in a uniform director, whereas the opposite occurs when the particle is in the zone of converging lines. The latter contribution is stronger than the former; hence, the overall free energy density change is negative; this means that the distortion entails an entropic gain, which translates into a negative value of the corresponding elastic constants combination, (K₁₁ − K₂₄).

FIG. 5.

View large Download slide

Contour plots showing the difference between the average excluded volume (a) in double twist or (b) in generalized double splay distortion with w = 2 and the corresponding quantity in a uniform director field. Calculations are performed for a pair of THSRs with M = 12 (ϕ = 0.282); one has its c.m. in the XZ plane and rotates around the other which has its c.m. at the origin of the reference frame $({\hat{n}}_{u} ‖ Z)$ ⁠. Color ranges from blue to red for increasingly negative values (white corresponds to zero). The deformations are defined in Table I, with qσ = 0.02. The gray lines in (b) are tangent to the splayed director field.

B. Elastic constants dependence on density and aspect ratio

Figures 6 and 7 show the elastic constants as a function of the packing fraction ϕ for THSRs with M = 6, 12, 24. Calculations were performed using Eqs. (20)–(22), (24), and (26), which are computationally advantageous over the calculation of deformation free energy density profiles described in Subsection IV A. The results obtained by the two methods coincide within numerical errors, as illustrated by the asterisks in the plots for the case M = 12. For K₁₃, we show also the values calculated according to Eq. (27), which are practically superimposed to those obtained from Eq. (26). For ease of visualization, the bend constants of all systems have been collected in the same plot, together with the corresponding S order parameters.

FIG. 6.

$Bend elastic constant as a function of the packing fraction ϕ, relative to the value at the isotropic–nematic coexistence, ϕN, for THSRs with M = 6 (dots), 12 (solid line) and 24 (dashed line). For M = 12, the value obtained from the analysis of the deformation free energy density reported in Fig. 3 is shown (asterisk). The inset shows the ϕ dependence of the nematic order parameter S for the three systems.$

View large Download slide

Bend elastic constant as a function of the packing fraction ϕ, relative to the value at the isotropic–nematic coexistence, ϕ_N, for THSRs with M = 6 (dots), 12 (solid line) and 24 (dashed line). For M = 12, the value obtained from the analysis of the deformation free energy density reported in Fig. 3 is shown (asterisk). The inset shows the ϕ dependence of the nematic order parameter S for the three systems.

FIG. 7.

$Elastic constants as a function of the packing fraction ϕ, for THSRs with (a) M = 6, (b) 12, and (c) 24. For K13, the solid line and the crosses represent values obtained from numerical evaluation of Eq. (26) and from Eq. (27), respectively. The dashed line corresponds to the renormalized splay elastic constant, defined in Eq. (32). In plot (b), the asterisks indicate the elastic coefficients obtained from the analysis of the deformation free energy density reported in Figs. 3 and 4.$

View large Download slide

Elastic constants as a function of the packing fraction ϕ, for THSRs with (a) M = 6, (b) 12, and (c) 24. For K₁₃, the solid line and the crosses represent values obtained from numerical evaluation of Eq. (26) and from Eq. (27), respectively. The dashed line corresponds to the renormalized splay elastic constant, defined in Eq. (32). In plot (b), the asterisks indicate the elastic coefficients obtained from the analysis of the deformation free energy density reported in Figs. 3 and 4.

As expected, with increasing particle length the isotropic–nematic transition occurs at lower ϕ and the bend constant increases dramatically, whereas the other two bulk constants, which are both smaller than K₃₃ and obey K₂₂ < K₁₁, are much less sensitive to the aspect ratio. What is more interesting, because much less known, is the behavior of the two surface-like constants: analogously to the bulk ones, they increase with increasing density and at any density they are higher than both K₁₁ and K₂₂, irrespective of the aspect ratio of the particles.

The results shown in Fig. 7 confirm the local instability of the uniform nematic state against double twist and certain variants of double splay as general features of hard rod nematics. This, however, does not necessarily translates into the appearance of a deformed ground state in finite samples because energetically favorable director distortions may be unable to fill 3D Euclidean space;⁵⁴ hence, they must be accompanied by defects or by other deformations, which give a positive contribution to the total free energy.⁵⁵ We must also mention that the value that we calculate for K₁₁ could be an underestimate, since an additional cost was pointed out for long linear particles, which originates from the coupling between changes in orientational order and density variations, and can be related to the decrease of entropy caused by the different density of particle heads and tails in the presence of splay.^56,57 Thus, a “renormalized” splay elastic constant can be approximated as the sum of the two contributions,

K_{11}^{R} \approx K_{11} + \frac{1}{2} k_{B} T l ρ_{0},

(32)

where l is the average longitudinal projection of a particle and ρ₀ = ρl is the areal number density of particles. The second term was found to dominate the splay constant in Monte Carlo simulations of relatively stiff polymers.⁵⁸ The dashed lines in the plots of Fig. 7 show the “renormalized” splay elastic constant of THSRs, calculated assuming

l = (M - 1) σ {⟨ | u_{Z} | ⟩}_{u}

⁠, and we can see that indeed this linear term gives a significant contribution to the splay constant, especially for long chains.

V. CONCLUSIONS

The main novelties of the present work can be summarized as follows:

We have presented a general formulation of Onsager theory for hard rod nematics that allows us to calculate the deformation free energy density for any arbitrary director field, without necessarily introducing the customary approximation of long wavelength distortions. The only assumption is a constant value of the nematic order parameter S, which might be relaxed in future, more general approaches.
From the analysis of Δ and double twist deformations, we obtain consistent estimates of the saddle–splay constant K₂₄. This result is intriguing because it implies an intrinsic tendency of hard rods towards double twist, i.e., towards spontaneous breaking of the mirror symmetry.⁴² For the same systems, we find also K₂₄ > K₁₁, but this result could be questionable, since the inequality would be reversed if the additional contribution to the splay constant coming from density-director coupling was taken into account.
The splay–bend contribution, which is missing in first order elastic theories, turns out to be crucial to explain the difference between the (microscopic) deformation free energy density calculated for variants of double splay, and the K₁₃ constant is found to be higher than all the other constants but the bend one.
We have obtained also explicit expressions for the full set of elastic constants and the results are in agreement with those obtained from direct analysis of the microscopic deformation free energy densities outlined in (i), in the limit of long wavelength deformations. Calculations for rods over a wide range of thermodynamic conditions confirm that K₂₄ > K₂₂ and large K₁₃ are general results, irrespective of aspect ratio and density.

The capability to obtain estimates of the full set of elastic constants is especially valuable for predicting and interpreting liquid crystal configurations in confined and structured geometries. In the future, it will be interesting to extend the approach, which here has been implemented for hard rods, to other forms of the interaction potential. In particular, it will be worth investigating whether the tendency towards spontaneous twist, associated with K₂₄ > K₂₂, is a distinct character of excluded volume interactions or a more general property, and how it is affected by features of the mesogenic particles, such as their shape and flexibility.

ACKNOWLEDGMENTS

We would like to thank Professor Kostas Daoulas, Professor Jonathan Selinger, Professor Tim Sluckin, and Professor Epifanio Virga for enlightening discussions and an anonymous reviewer for insightful suggestions. We acknowledge the CAPRI Initiative (Calcolo ad Alte Prestazioni per la Ricerca e l’Innovazione, University of Padova Strategic Research Infrastructure Grant No. 2017) for HPC resources and technical support, and the contribution of the COST Action Grant No. CA17139. D.R. gratefully acknowledges Fondazione CARIPARO for funding his Ph.D. scholarship.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Davide Revignas: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal). Alberta Ferrarini: Conceptualization (equal); Funding acquisition (lead); Methodology (equal); Resources (equal); Supervision (lead); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A: MACROSCOPIC-TO-MICROSCOPIC CORRESPONDENCE IN THE FREE ENERGY OF DEFORMATION

In Ref. 42, it was shown how to link the integral of the elastic free energy density [Eq. ] performed over a sub-region of volume V contained in a nematic sample of volume

V

to the following microscopic expression for the deformation free energy:

\begin{align} \int_{V} d R Δ a_{d}^{def} (R, q) = - \frac{η ρ^{2} k_{B} T}{2} \int_{V} d R_{A} \int_{V} d R_{B} \int d {\hat{u}}_{A} \int d {\hat{u}}_{B} \\ \times f_{d} ({\hat{u}}_{A}, R_{A}, q) f_{d} ({\hat{u}}_{B}, R_{B}, q) \\ \times e_{A B} (R_{A B}, {\hat{u}}_{A}, {\hat{u}}_{B}) - (V a_{u}^{ex}) . \end{align}

(A1)

The first term on the rhs of Eq. comes from the division of the original integration domain

V

for the variable R_A into smaller sub-regions separated by virtual walls. In particular, the sub-region of volume V we are interested in is thought of as being far enough from the nematic sample’s physical boundaries that its free energy is purely bulk in nature. The second term at the rhs of Eq. is simply the free energy of the region when the nematic director field is uniform.

Equation establishes an unambiguous macroscopic (first member) to microscopic (second member) connection. In principle, there is no lower bound to the size of the volume of interest V: it is therefore tempting to push the macroscopic-to-microscopic connection even further by shrinking such domain down to a single point. With this operation, a correspondence between R and R_A is found, and by dropping the now superfluous A subscript one can write

\begin{align} Δ a_{d}^{def} (R, q) = - \frac{η ρ^{2} k_{B} T}{2} \times \int d {\hat{u}}_{A} \int d {\hat{u}}_{B} f_{d} ({\hat{u}}_{A}, R, q) f_{d} ({\hat{u}}_{B}, R_{B}, q) \\ \times \int_{V} d R_{B} e_{A B} (R_{A B}, {\hat{u}}_{A}, {\hat{u}}_{B}) - a_{u}^{ex} . \end{align}

(A2)

Equation is a more general form of Eqs. –, valid for any point R.

In the past, an objection was raised³² on the legitimacy of this mapping to a density: indeed, it was argued that, by performing changes of variables involving both R_A and R_B in Eq. , one could obtain different expressions for

Δ a_{d}^{def}

that in turn would lead to different values for the surface-like elastic constants. As an example, the following change of variables was proposed:

\{\begin{cases} R_{0} = (R_{A} + R_{B}) / 2 \\ R_{A B} = R_{B} - R_{A} . \end{cases}

(A3)

By applying this transformation in the first term at the rhs of Eq. , both in the kernels and in the integration domains, one can in principle establish this new macroscopic-to-microscopic correspondence

\begin{align} \int_{\tilde{V}} d R Δ {\tilde{a}}_{d}^{def} (R, q) = - \frac{η ρ^{2} k_{B} T}{2} \int_{\tilde{V}} d R_{0} \int_{\tilde{V} (R_{0})} d R_{A B} \\ \times \int d {\hat{u}}_{A} \int d {\hat{u}}_{B} f_{d} ({\hat{u}}_{A}, R_{0} - R_{A B} / 2, q) \\ \times f_{d} ({\hat{u}}_{B}, R_{0} + R_{A B} / 2, q) e_{A B} (R_{A B}, {\hat{u}}_{A}, {\hat{u}}_{B}) \\ - (\tilde{V} a_{u}^{ex}) . \end{align}

(A4)

Now, by shrinking again the transformed volume

\tilde{V}

down to a point, one would obtain a different expression for the free energy density,

Δ {\tilde{a}}_{d}^{def}

⁠. However,

Δ {\tilde{a}}_{d}^{def}

has the unphysical characteristic of being R₀-dependent even in absence of deformation, through its domain of integration. Analogous conclusions were derived, from a different perspective, in Ref. 33.⁵⁹

APPENDIX B: EQUIVALENCES BETWEEN INTEGRALS

In this Appendix, we derive some equivalences between integrals employed in the main text, which are a consequence of the uniaxial symmetry of the undeformed ODF [Eq. ]. There are two groups of such equivalences: the first one is

⟨E_{R_{X}^{2}} u_{B X}^{2} u_{B Z}^{2}⟩ = ⟨E_{R_{Y}^{2}} u_{B Y}^{2} u_{B Z}^{2}⟩,

(B1)

⟨E_{R_{X}^{2}} (u_{B X}^{2} - u_{B Z}^{2})⟩ = ⟨E_{R_{Y}^{2}} (u_{B Y}^{2} - u_{B Z}^{2})⟩,

(B2)

⟨E_{R_{X} R_{Z}} u_{B X} u_{B Z}⟩ = ⟨E_{R_{Y} R_{Z}} u_{B Y} u_{B Z}⟩,

(B3)

⟨E_{R_{X}^{2}} u_{B Z}^{2}⟩ = ⟨E_{R_{Y}^{2}} u_{B Z}^{2}⟩ .

(B4)

The second group is

2 ⟨E_{R_{X} R_{Y}} u_{B X} u_{B Y}⟩ = ⟨E_{R_{X}^{2}} u_{B X}^{2}⟩ - ⟨E_{R_{Y}^{2}} u_{B X}^{2}⟩,

(B5)

2 ⟨E_{R_{X} R_{Y}} u_{B X} u_{B Y} u_{B Z}^{2}⟩ = ⟨E_{R_{X}^{2}} u_{B X}^{2} u_{B Z}^{2}⟩ - ⟨E_{R_{Y}^{2}} u_{B X}^{2} u_{B Z}^{2}⟩ .

(B6)

We will explicitly demonstrate one equivalence per group since the procedure is similar for all the others.

Let us start from the lhs of Eq. , whose explicit form is

\begin{align} ⟨E_{R_{X}^{2}} u_{B X}^{2} u_{B Z}^{2}⟩ = \frac{1}{Q^{2}} \int d {\hat{u}}_{A} e^{c {[{\hat{u}}_{A} \cdot \hat{Z}]}^{2}} \int d {\hat{u}}_{B} e^{c {[{\hat{u}}_{B} \cdot \hat{Z}]}^{2}} \\ \times {[{\hat{u}}_{B} \cdot \hat{X}]}^{2} {[{\hat{u}}_{B} \cdot \hat{Z}]}^{2} \int d R e_{A B} \\ \times ({\hat{u}}_{A}, {\hat{u}}_{B}, R) {[R \cdot \hat{X}]}^{2} . \end{align}

(B7)

If we perform a rotation of 90° around the Z axis of the reference frame, the quantities involved in the last expression transform according to the following relations:

{\hat{u}}_{A} \to {\hat{u}}_{A}^{'} \hat{X} \to \hat{Y},

(B8)

{\hat{u}}_{B} \to {\hat{u}}_{B}^{'} \hat{Y} \to - \hat{X},

(B9)

R_{A B} \to R_{A B}^{'} \hat{Z} \to \hat{Z} .

(B10)

Therefore, we obtain

\begin{align} \frac{1}{Q^{2}} \int d {\hat{u}}_{A}^{'} e^{c {[{\hat{u}}_{A}^{'} \cdot \hat{Z}]}^{2}} \int d {\hat{u}}_{B}^{'} e^{c {[{\hat{u}}_{B}^{'} \cdot \hat{Z}]}^{2}} {[- {\hat{u}}_{B}^{'} \cdot \hat{Y}]}^{2} {[{\hat{u}}_{B}^{'} \cdot \hat{Z}]}^{2} \\ \times \int d R^{'} e_{A B} ({\hat{u}}_{A}^{'}, {\hat{u}}_{B}^{'}, R^{'}) {[- R^{'} \cdot \hat{Y}]}^{2} = ⟨E_{{R^{'}}_{Y}^{2}} {u^{'}}_{B Y}^{2} {u^{'}}_{B Z}^{2}⟩ . \end{align}

(B11)

By dropping the prime, one recovers the rhs of Eq. .

Let us now examine the lhs of Eq. from the second group of equivalences

\begin{align} 2 ⟨E_{R_{X} R_{Y}} u_{B X} u_{B Y}⟩ = \frac{2}{Q^{2}} \int d {\hat{u}}_{A} e^{c {[{\hat{u}}_{A} \cdot \hat{Z}]}^{2}} \int d {\hat{u}}_{B} e^{c {[{\hat{u}}_{B} \cdot \hat{Z}]}^{2}} \\ \times [{\hat{u}}_{B} \cdot \hat{X}] [{\hat{u}}_{B} \cdot \hat{Y}] \int d R e_{A B} ({\hat{u}}_{A}, {\hat{u}}_{B}, R) \\ \times [R \cdot \hat{X}] [R \cdot \hat{Y}] . \end{align}

(B12)

If we perform a rotation of 45° around the Z axis of the reference frame, the quantities involved in this expression transform according to the following relations:

{\hat{u}}_{A} \to {\hat{u}}_{A}^{'} \hat{X} \to \frac{1}{\sqrt{2}} (\hat{X} - \hat{Y}),

(B13)

{\hat{u}}_{B} \to {\hat{u}}_{B}^{'} \hat{Y} \to \frac{1}{\sqrt{2}} (\hat{X} + \hat{Y}),

(B14)

R_{A B} \to R_{A B}^{'} \hat{Z} \to \hat{Z} .

(B15)

Therefore, Eq. becomes

\begin{align} 2 ⟨E_{R_{X} R_{Y}} u_{B X} u_{B Y}⟩ = \frac{1}{2 Q^{2}} \int d {\hat{u}}_{A}^{'} e^{c {[{\hat{u}}_{A}^{'} \cdot \hat{Z}]}^{2}} \int d {\hat{u}}_{B}^{'} e^{c {[{\hat{u}}_{B}^{'} \cdot \hat{Z}]}^{2}} \\ \times [{\hat{u}}_{B}^{'} \cdot (\hat{X} - \hat{Y})] [{\hat{u}}_{B}^{'} \cdot (\hat{X} + \hat{Y})] \\ \times \int d R^{'} e_{A B} ({\hat{u}}_{A}^{'}, {\hat{u}}_{B}^{'}, R^{'}) [R^{'} \cdot (\hat{X} - \hat{Y})] \\ \times [R^{'} \cdot (\hat{X} + \hat{Y})] . \end{align}

(B16)

Rearranging the terms, we obtain

\begin{align} 4 ⟨E_{R_{X} R_{Y}} u_{B X} u_{B Y}⟩ = \frac{1}{Q^{2}} \int d {\hat{u}}_{A}^{'} e^{c {[{\hat{u}}_{A}^{'} \cdot \hat{Z}]}^{2}} \int d {\hat{u}}_{B}^{'} e^{c {[{\hat{u}}_{B}^{'} \cdot \hat{Z}]}^{2}} \\ \times \int d R^{'} e_{A B} ({\hat{u}}_{A}^{'}, {\hat{u}}_{B}^{'}, R^{'}) \\ \times [u_{B X}^{'} - u_{B Y}^{'}] [u_{B X}^{'} + u_{B Y}^{'}] [R_{B X}^{'} - R_{B Y}^{'}] [R_{B X}^{'} + R_{B Y}^{'}] \\ = \frac{1}{Q^{2}} \int d {\hat{u}}_{A}^{'} e^{c {[{\hat{u}}_{A}^{'} \cdot \hat{Z}]}^{2}} \int d {\hat{u}}_{B}^{'} e^{c {[{\hat{u}}_{B}^{'} \cdot \hat{Z}]}^{2}} \\ \times \int d R^{'} e_{A B} ({\hat{u}}_{A}^{'}, {\hat{u}}_{B}^{'}, R^{'}) [{u^{'}}_{B X}^{2} - {u^{'}}_{B Y}^{2}] \\ \times [{R^{'}}_{B X}^{2} - {R^{'}}_{B Y}^{2}], \end{align}

(B17)

and then by dropping the prime, we can write

\begin{align} 4 ⟨E_{R_{X} R_{Y}} u_{B X} u_{B Y}⟩ & = ⟨E_{R_{X}^{2}} u_{B X}^{2}⟩ + ⟨E_{R_{Y}^{2}} u_{B Y}^{2}⟩ + \\ - ⟨E_{R_{X}^{2}} u_{B Y}^{2}⟩ - ⟨E_{R_{Y}^{2}} u_{B X}^{2}⟩ . \end{align}

(B18)

Using the equivalences

⟨E_{R_{Y}^{2}} u_{B Y}^{2}⟩ = ⟨E_{R_{X}^{2}} u_{B X}^{2}⟩

and

⟨E_{R_{Y}^{2}} u_{B X}^{2}⟩ = ⟨E_{R_{X}^{2}} u_{B Y}^{2}⟩

⁠, which again can be obtained performing 90° around the Z axis of the reference frame, Eq. is recovered.

REFERENCES

1.

C. W.

Oseen

, “

The theory of liquid crystals

,”

Trans. Faraday Soc.

29

,

883

–

899

(

1933

).

https://doi.org/10.1039/tf9332900883

Google Scholar

Crossref

2.

F. C.

Frank

, “

I. Liquid crystals. On the theory of liquid crystals

,”

Discuss. Faraday Soc.

25

,

19

–

28

(

1958

).

https://doi.org/10.1039/df9582500019

Google Scholar

Crossref

3.

D. W.

Allender

,

G. P.

Crawford

, and

J. W.

Doane

, “

Determination of the liquid-crystal surface elastic constant K₂₄

,”

Phys. Rev. Lett.

67

,

1442

–

1445

(

1991

).

https://doi.org/10.1103/physrevlett.67.1442

Google Scholar

Crossref

PubMed

4.

O. D.

Lavrentovich

and

V. M.

Pergamenshchik

, “

Stripe domain phase of a thin nematic film and the K₁₃ divergence term

,”

Phys. Rev. Lett.

73

,

979

–

982

(

1994

).

https://doi.org/10.1103/physrevlett.73.979

Google Scholar

Crossref

PubMed

5.

R. D.

Polak

,

G. P.

Crawford

,

B. C.

Kostival

,

J. W.

Doane

, and

S.

Žumer

, “

Optical determination of the saddle-splay elastic constant K₂₄ in nematic liquid crystals

,”

Phys. Rev. B

49

,

R978

–

R981

(

1994

).

https://doi.org/10.1103/physreve.49.r978

Google Scholar

Crossref

6.

J.

Nehring

and

A.

Saupe

, “

On the elastic theory of uniaxial liquid crystals

,”

J. Chem. Phys.

54

,

337

–

343

(

1971

).

https://doi.org/10.1063/1.1674612

Google Scholar

Crossref

7.

J. L.

Ericksen

, “

Inequalities in liquid crystal theory

,”

Phys. Fluids

9

,

1205

–

1207

(

1966

).

https://doi.org/10.1063/1.1761821

Google Scholar

Crossref

8.

Ž.

Kos

and

M.

Ravnik

, “

Relevance of saddle-splay elasticity in complex nematic geometries

,”

Soft Matter

12

,

1313

–

1323

(

2016

).

https://doi.org/10.1039/c5sm02417j

Google Scholar

Crossref

PubMed

9.

H.

Kikuchi

, “

Liquid crystalline blue phases

,” in

Structure and Bonding - Liquid Crystalline Functional Assemblies and Their Supramolecular Structures

, edited by

T.

Kato

(

Springer

,

Berlin

,

2007

), pp.

99

–

117

.

Google Scholar

Crossref

10.

J. V.

Selinger

, “

Interpretation of saddle-splay and the Oseen-Frank free energy in liquid crystals

,”

Liq. Cryst. Rev.

6

,

129

–

142

(

2018

).

https://doi.org/10.1080/21680396.2019.1581103

Google Scholar

Crossref

11.

T.

Machon

and

G. P.

Alexander

, “

Umbilic lines in orientational order

,”

Phys. Rev. X

6

,

011033

(

2016

).

https://doi.org/10.1103/physrevx.6.011033

Google Scholar

Crossref

12.

G.

Barbero

and

L. R.

Evangelista

,

An Elementary Course on the Continuum Theory for Nematic Liquid Crystals

(

World Scientific

,

Singapore

,

2001

).

Google Scholar

13.

H.

Zocher

, “

The effect of a magnetic field on the nematic state

,”

Trans. Faraday Soc.

29

,

945

–

957

(

1933

).

https://doi.org/10.1039/tf9332900945

Google Scholar

Crossref

14.

V. M.

Pergamenshchik

, “

Phenomenological approach to the problem of the K₁₃ surfacelike elastic term in the free energy of a nematic liquid crystal

,”

Phys. Rev. E

48

,

1254

–

1264

(

1993

).

https://doi.org/10.1103/physreve.48.1254

Google Scholar

Crossref

15.

J. P.

De Gennes

and

J.

Prost

,

The Physics of Liquid Crystals

(

Oxford University Press

,

Oxford

,

1993

).

Google Scholar

Crossref

16.

N. V.

Madhusudana

and

R.

Pratibha

, “

Experimental determination of the elastic constant K₁₃ of a nematic liquid crystal

,”

Mol. Cryst. Liq. Cryst.

179

,

207

–

216

(

1990

).

https://doi.org/10.1080/00268949008055371

Google Scholar

Crossref

17.

A.

Sparavigna

,

O. D.

Lavrentovich

, and

A.

Strigazzi

, “

Periodic stripe domains and hybrid-alignment regime in nematic liquid crystals: Threshold analysis

,”

Phys. Rev. E

49

,

1344

–

1352

(

1994

).

https://doi.org/10.1103/physreve.49.1344

Google Scholar

Crossref

18.

G. P.

Crawford

,

D. W.

Allender

,

J. W.

Doane

,

M.

Vilfan

, and

I.

Vilfan

, “

Finite molecular anchoring in the escaped-radial nematic configuration: A H2-NMR study

,”

Phys. Rev. A

44

,

2570

–

2577

(

1991

).

https://doi.org/10.1103/physreva.44.2570

Google Scholar

Crossref

PubMed

19.

G. P.

Crawford

,

D. W.

Allender

, and

J. W.

Doane

, “

Surface elastic and molecular-anchoring properties of nematic liquid crystals confined to cylindrical cavities

,”

Phys. Rev. A

45

,

8693

–

8708

(

1992

).

https://doi.org/10.1103/physreva.45.8693

Google Scholar

Crossref

PubMed

20.

R. J.

Ondris-Crawford

,

G. P.

Crawford

,

S.

Zumer

, and

J. W.

Doane

, “

Curvature-induced configuration transition in confined nematic liquid crystals

,”

Phys. Rev. Lett.

70

,

194

–

197

(

1993

).

https://doi.org/10.1103/physrevlett.70.194

Google Scholar

Crossref

PubMed

21.

Y.

Xia

,

A. A.

DeBenedictis

,

D. S.

Kim

,

S.

Chen

,

S.-U.

Kim

,

D. J.

Cleaver

,

T. J.

Atherton

, and

S.

Yang

, “

Programming emergent symmetries with saddle-splay elasticity

,”

Nat. Commun.

10

,

5104

(

2019

).

https://doi.org/10.1038/s41467-019-13012-9

Google Scholar

Crossref

PubMed

22.

E.

Pairam

,

J.

Vallamkondu

,

V.

Koning

,

B. C.

van Zuiden

,

P. W.

Ellis

,

M. A.

Bates

,

V.

Vitelli

, and

A.

Fernandez-Nieves

, “

Stable nematic droplets with handles

,”

Proc. Natl. Acad. Sci. U. S. A.

110

,

9295

–

9300

(

2013

).

https://doi.org/10.1073/pnas.1221380110

Google Scholar

Crossref

PubMed

23.

D. S.

Miller

and

N. L.

Abbott

, “

Influence of droplet size, ph and ionic strength on endotoxin-triggered ordering transitions in liquid crystalline droplets

,”

Soft Matter

9

,

374

–

382

(

2013

).

https://doi.org/10.1039/c2sm26811f

Google Scholar

Crossref

PubMed

24.

J.

Lydon

, “

Chromonic review

,”

J. Mater. Chem.

20

,

10071

–

10099

(

2010

).

https://doi.org/10.1039/b926374h

Google Scholar

Crossref

25.

Z. S.

Davidson

,

L.

Kang

,

J.

Jeong

,

T.

Still

,

P. J.

Collings

,

T. C.

Lubensky

, and

A. G.

Yodh

, “

Chiral structures and defects of lyotropic chromonic liquid crystals induced by saddle-splay elasticity

,”

Phys. Rev. E

91

,

050501

(

2015

).

https://doi.org/10.1103/physreve.91.050501

Google Scholar

Crossref

26.

K.

Nayani

,

R.

Chang

,

J.

Fu

,

P. W.

Ellis

,

A.

Fernandez-Nieves

,

J. O.

Park

, and

M.

Srinivasarao

, “

Spontaneous emergence of chirality in achiral lyotropic chromonic liquid crystals confined to cylinders

,”

Nat. Commun.

6

,

8067

(

2015

).

https://doi.org/10.1038/ncomms9067

Google Scholar

Crossref

PubMed

27.

J.

Eun

,

S.-J.

Kim

, and

J.

Jeong

, “

Effects of chiral dopants on double-twist configurations of lyotropic chromonic liquid crystals in a cylindrical cavity

,”

Phys. Rev. E

100

,

012702

(

2019

).

https://doi.org/10.1103/physreve.100.012702

Google Scholar

Crossref

PubMed

28.

A.

Javadi

,

J.

Eun

, and

J.

Jeong

, “

Cylindrical nematic liquid crystal shell: Effect of saddle-splay elasticity

,”

Soft Matter

14

,

9005

–

9011

(

2018

).

https://doi.org/10.1039/c8sm01829d

Google Scholar

Crossref

PubMed

29.

T. C.

Lubensky

, “

Confined chromonics and viral membranes

,”

Mol. Cryst. Liq. Cryst.

646

,

235

–

241

(

2017

).

https://doi.org/10.1080/15421406.2017.1288010

Google Scholar

Crossref

30.

S.

Fumeron

,

F.

Moraes

, and

E.

Pereira

, “

Retrieving the saddle-splay elastic constant K₂₄ of nematic liquid crystals from an algebraic approach

,”

Eur. Phys. J. E

39

,

83

(

2016

).

https://doi.org/10.1140/epje/i2016-16083-8

Google Scholar

Crossref

31.

J.

Nehring

and

A.

Saupe

, “

Calculation of the elastic constants of nematic liquid crystals

,”

J. Chem. Phys.

56

,

5527

–

5528

(

1972

).

https://doi.org/10.1063/1.1677071

Google Scholar

Crossref

32.

A. M.

Somoza

and

P.

Tarazona

, “

Density functional theory of the elastic constants of a nematic liquid crystal

,”

Mol. Phys.

72

,

911

–

926

(

1991

).

https://doi.org/10.1080/00268979100100641

Google Scholar

Crossref

33.

P. I. C.

Teixeira

,

V. M.

Pergamenshchik

, and

T. J.

Sluckin

, “

A model calculation of the surface elastic constants of a nematic liquid crystal

,”

Mol. Phys.

80

,

1339

–

1357

(

1993

).

https://doi.org/10.1080/00268979300103071

Google Scholar

Crossref

34.

G.

Barbero

and

C.

Oldano

, “

Distortion in the surface layers of nematic liquid crystals: An elastic constant approach

,”

Mol. Cryst. Liq. Cryst.

170

,

99

–

106

(

1989

).

https://doi.org/10.1080/00268948908047751

Google Scholar

Crossref

35.

J.

Stelzer

,

L.

Longa

, and

H. R.

Trebin

, “

Molecular dynamics simulations of a Gay–Berne nematic liquid crystal: Elastic properties from direct correlation functions

,”

J. Chem. Phys.

103

,

3098

–

3107

(

1995

).

https://doi.org/10.1063/1.470268

Google Scholar

Crossref

36.

M.

Cestari

and

A.

Ferrarini

, “

Curvature elasticity of nematic liquid crystals: Simply a matter of molecular shape? Insights from atomistic modeling

,”

Soft Matter

5

,

3879

–

3887

(

2009

).

https://doi.org/10.1039/b904606b

Google Scholar

Crossref

37.

A. A.

Joshi

,

J. K.

Whitmer

,

O.

Guzmán

,

N. L.

Abbott

, and

J. J.

de Pablo

, “

Measuring liquid crystal elastic constants with free energy perturbations

,”

Soft Matter

10

,

882

–

893

(

2014

).

https://doi.org/10.1039/c3sm51919h

Google Scholar

Crossref

PubMed

38.

H.

Sidky

,

J. J.

de Pablo

, and

J. K.

Whitmer

, “

In silico measurement of elastic moduli of nematic liquid crystals

,”

Phys. Rev. Lett.

120

,

107801

(

2018

).

https://doi.org/10.1103/physrevlett.120.107801

Google Scholar

Crossref

PubMed

39.

M.

Cestari

,

A.

Bosco

, and

A.

Ferrarini

, “

Molecular field theory with atomistic modeling for the curvature elasticity of nematic liquid crystals

,”

J. Chem. Phys.

131

,

054104

(

2009

).

https://doi.org/10.1063/1.3193555

Google Scholar

Crossref

PubMed

40.

M. A.

Osipov

and

G.

Pajak

, “

Effect of polar intermolecular interactions on the elastic constants of bent-core nematics and the origin of the twist-bend phase

,”

Eur. Phys. J. E

39

,

45

(

2016

).

https://doi.org/10.1140/epje/i2016-16045-2

Google Scholar

Crossref

41.

D.

Revignas

and

A.

Ferrarini

, “

From bend to splay dominated elasticity in nematics

,”

Crystals

11

,

831

(

2021

).

https://doi.org/10.3390/cryst11070831

Google Scholar

Crossref

42.

D.

Revignas

and

A.

Ferrarini

, “

Spontaneous twisting of achiral hard rod nematics

,”

Phys. Rev. Lett.

130

,

028102

(

2023

).

https://doi.org/10.1103/PhysRevLett.130.028102

Google Scholar

Crossref

PubMed

43.

L.

Onsager

, “

The effects of shape on the interaction of colloidal particles

,”

Ann. N. Y. Acad. Sci.

51

,

627

–

659

(

1949

).

https://doi.org/10.1111/j.1749-6632.1949.tb27296.x

Google Scholar

Crossref

44.

J. P.

Straley

, “

Frank elastic constants of the hard-rod liquid crystal

,”

Phys. Rev. A

8

,

2181

–

2183

(

1973

).

https://doi.org/10.1103/physreva.8.2181

Google Scholar

Crossref

45.

S.

Varga

and

I.

Szalai

, “

Modified Parsons-Lee theory for fluids of linear fused hard sphere chains

,”

Mol. Phys.

98

,

693

–

698

(

2000

).

https://doi.org/10.1080/00268970009483337

Google Scholar

Crossref

46.

J. D.

Parsons

, “

Nematic ordering in a system of rods

,”

Phys. Rev. A

19

,

1225

–

1230

(

1979

).

https://doi.org/10.1103/physreva.19.1225

Google Scholar

Crossref

47.

S.-D.

Lee

, “

A numerical investigation of nematic ordering based on a simple hard-rod model

,”

J. Chem. Phys.

87

,

4972

–

4974

(

1987

).

https://doi.org/10.1063/1.452811

Google Scholar

Crossref

48.

D. A.

McQuarrie

,

Statistical Mechanics

(

Harper & Row

,

New York

,

1976

).

Google Scholar

49.

S.

Paparini

and

E.

Virga

, “

An elastic quartic twist theory for chromonic liquid crystals

,”

J. Elasticity

(published online, 2023).

https://doi.org/10.1007/s10659-022-09983-4

50.

M. F.

Sanner

,

A. J.

Olson

, and

J.-C.

Spehner

, “

Reduced surface: An efficient way to compute molecular surfaces

,”

Biopolymers

38

,

305

–

320

(

1996

).

https://doi.org/10.1002/(sici)1097-0282(199603)38:3<305::aid-bip4>3.0.co;2-y

Google Scholar

Crossref

PubMed

51.

W. H.

Press

,

B. P.

Flannery

,

S. A.

Teukolsky

, and

W. T.

Vetterling

,

Numerical Recipes: The Art of Scientific Computing

(

Cambridge University Press

,

Cambridge

,

1986

).

Google Scholar

52.

J. C.

Lagarias

,

J. A.

Reeds

,

M. H.

Wright

, and

P. E.

Wright

, “

Convergence properties of the Nelder–Mead simplex method in low dimensions

,”

SIAM J. Optim.

9

,

112

–

147

(

1998

).

https://doi.org/10.1137/s1052623496303470

Google Scholar

Crossref

53.

MATLAB, 9.2.0.538062 (R2017a),

Natick

,

MA

,

2017

.

54.

J. V.

Selinger

, “

Director deformations, geometric frustration, and modulated phases in liquid crystals

,”

Annu. Rev. Condens. Matter Phys.

49

,

71

(

2022

).

https://doi.org/10.1146/annurev-conmatphys-031620-105712

Google Scholar

Crossref

55.

J. P.

Sethna

,

D. C.

Wright

, and

N. D.

Mermin

, “

Relieving cholesteric frustration: The blue phase in a curved space

,”

Phys. Rev. Lett.

51

,

467

–

470

(

1983

).

https://doi.org/10.1103/physrevlett.51.467

Google Scholar

Crossref

56.

R. B.

Meyer

, in

Polymer Liquid Crystals

, edited by

A.

Ciferri

,

W. R.

Kringbaum

, and

R. B.

Meyer

(

Academic Press

,

New York

,

1982

).

Google Scholar

57.

P.

Le Doussal

and

D. R.

Nelson

, “

Statistical mechanics of directed polymer melts

,”

Europhys. Lett.

15

,

161

–

166

(

1991

).

https://doi.org/10.1209/0295-5075/15/2/009

Google Scholar

Crossref

58.

P.

Gemünden

and

K. C.

Daoulas

, “

Fluctuation spectra in polymer nematics and Frank elastic constants: A coarse-grained modelling study

,”

Soft Matter

11

,

532

–

544

(

2015

).

https://doi.org/10.1039/c4sm02075h

Google Scholar

Crossref

PubMed

59.

Changes of variables of the form of Eq. (A3) were sometimes employed to calculate the bulk elastic constants, and for this case we numerically verified the invariance of the results.

Published open access through an agreement with Universita degli Studi di Padova Dipartimento di Scienze Chimiche

On the elusive saddle–splay and splay–bend elastic constants of nematic liquid crystals

I. INTRODUCTION

II. THEORY

A. Macroscopic deformation analysis

B. Microscopic free energy density of deformation

C. Explicit expressions for the elastic constants

III. NUMERICAL IMPLEMENTATION

A. Model details

B. Calculation of integrals

C. Isotropic–nematic coexistence and equilibrium order parameter

D. Calculations for the contour plots of the average excluded volume

IV. RESULTS

A. Microscopic free energy density for different deformations

B. Elastic constants dependence on density and aspect ratio

V. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: MACROSCOPIC-TO-MICROSCOPIC CORRESPONDENCE IN THE FREE ENERGY OF DEFORMATION

APPENDIX B: EQUIVALENCES BETWEEN INTEGRALS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

On the elusive saddle–splay and splay–bend elastic constants of nematic liquid crystals

I. INTRODUCTION

II. THEORY

A. Macroscopic deformation analysis

B. Microscopic free energy density of deformation

C. Explicit expressions for the elastic constants

III. NUMERICAL IMPLEMENTATION

A. Model details

B. Calculation of integrals

C. Isotropic–nematic coexistence and equilibrium order parameter

D. Calculations for the contour plots of the average excluded volume

IV. RESULTS

A. Microscopic free energy density for different deformations

B. Elastic constants dependence on density and aspect ratio

V. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: MACROSCOPIC-TO-MICROSCOPIC CORRESPONDENCE IN THE FREE ENERGY OF DEFORMATION

APPENDIX B: EQUIVALENCES BETWEEN INTEGRALS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only