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*_{24}) and splay–bend (*K*_{13}) 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*_{24} is larger than both *K*_{11} and *K*_{22} over a wide range of particle lengths and densities. Moreover, the *K*_{13} 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

**) dependent director, $n\u0302(R)$,**

*R*^{1,2}

*K*

_{11},

*K*

_{22}, and

*K*

_{33}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. (1).

^{3–5}In other cases, the following alternative form has been employed:

^{6,7}

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.^{8} 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*_{24} < min(*K*_{11}, *K*_{22}).^{7} In particular, *K*_{22} − *K*_{24} < 0 would promote a double twist deformation, as in blue phases.^{9}

In a recent reformulation of nematic elasticity,^{10} 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*_{24}. Thus, a form of the free energy density is proposed, with the saddle–splay contribution explicitly regarded as a bulk term.

^{10,12}It was neglected by Frank,

^{2}but it was included by Oseen,

^{1}Zocher,

^{13}and Nehring and Saupe.

^{6}This term is generally denoted as (

*mixed*)

*splay*–

*bend*, and its contribution to the free energy density can be expressed as

*K*

_{13}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.

^{14}

Bulk elastic constants can be measured by independently exciting the relative modes, through a suitable use of sample geometry, boundary conditions and external fields.^{15} 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*_{13}, comparable with the bulk constants, was determined in hybrid nematic cells,^{16} whereas for 5CB *K*_{13} ≈ −0.20*K*_{11} was inferred from analysis of the periodicity of the strip domain phase of a hybrid aligned nematic film with azimuthally degenerate boundary conditions.^{4} For the saddle–splay constant *K*_{24} some more data, yet rather controversial, are available. *K*_{24} 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^{5} 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*_{24} was proposed, based on the use of combined topographic and chemical surface patterning to create saddle–splay distortions: *K*_{24} ≈ 0.5*K*_{22} was estimated for 5CB and for 8CB.^{21} 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*_{24} ≈ *K*_{22}^{22} was inferred from the formation of doubly twisted configurations in toroidal droplets, whereas *K*_{24} ≈ 0.7*K*_{11} was estimated from the analysis of the radial to bipolar ordering transition, induced by endotoxin in spherical droplets.^{23} 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.^{24} 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,^{9} when confined in cylindrical capillaries with degenerate planar boundary conditions^{25–27} and in cylindrical shells.^{28} 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*_{24} values were inferred from the analysis of the director configurations in these systems: *K*_{24}/*K*_{22} = 27.5^{25} and 3.25^{29,30} for SSY, *K*_{24}/*K*_{22} = 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),^{31} and the relationships *K*_{11}:*K*_{22}:*K*_{33}:*K*_{13} = 5:11:5:−6 were obtained, together with very small positive *K*_{24}. Remarkably, the very possibility of unambiguously defining *K*_{13} and *K*_{24} from a microscopic theory was questioned.^{32} This objection was discussed in a subsequent study,^{33} where the surface-like elastic constants of a Gay–Berne fluid were evaluated using density functional theory. *K*_{13} and *K*_{24} close to *K*_{22}/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*_{33} = *K*_{11}. 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*_{24} close to *K*_{22} and small negative *K*_{13} 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*_{24} > *K*_{22} for PAA, and *K*_{24} ≈ *K*_{22} for 5CB and 8CB, whereas *K*_{13} 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*_{24} ≈ *K*_{22} was obtained,^{37} whereas for an atomistic model of 5CB, *K*_{24} between 1.5 and 2.5 pN, smaller than *K*_{22}, was found.^{38}

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*_{22} < *K*_{11} < *K*_{33}, but also in less conventional systems.^{39–41} Recently,^{42} 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*_{24} 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*_{24} > *K*_{22} 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*_{13}. 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

**(expressed with respect to a fixed laboratory frame**

*R**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 $n\u0302u\Vert Z\u0302$ is recovered; and (iii) the director at the origin is in any case parallel to the laboratory

*Z*axis, $n\u0302d(0,q)=Z\u0302$. Using these expressions, the elastic free energy density, sum of Eqs. (1) and (3), at the origin of the reference frame,

**= 0, has a parabolic dependence on**

*R**q*and the parabolic coefficients are given by linear combinations, $Cd$, of elastic constants

*K*

_{24}, whereas the other two, double twist and double splay [see Fig. 1(b)], are related to linear combinations of

*K*

_{24},

*K*

_{11}, and

*K*

_{22}. 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*

_{13}constant.

Deformation name (d) . | $n\u0302d(R,q)$ . | $Cd=\Delta addef(0,q)/q2$ . |
---|---|---|

Splay (S) | $[qRX,0,1+qRZ]/q2RX2+(1+qRZ)2$ | K_{11}/2 |

Twist (T) | [sin(qR_{Y}), 0, cos(qR_{Y})] | K_{22}/2 |

Bend (B) | $[qRZ,0,1\u2212qRX]/q2RZ2+(1\u2212qRX)2$ | K_{33}/2 |

Δ deformation (Δ) | $[qRX,\u2212qRY,1]/q2(RX2+RY2)+1$ | 2K_{24} |

Double twist (DT) | $[\u2212qRY,qRX,1]/q2(RX2+RY2)+1$ | 2(K_{22} − K_{24}) |

Generalized | $[qRX,qRY,1+qwRZ]/q2(RX2+RY2)+(1+qwRZ)2$ | 2[K_{11} + (2 − w)K_{13} − K_{24}]a |

double splay (DS) |

Deformation name (d) . | $n\u0302d(R,q)$ . | $Cd=\Delta addef(0,q)/q2$ . |
---|---|---|

Splay (S) | $[qRX,0,1+qRZ]/q2RX2+(1+qRZ)2$ | K_{11}/2 |

Twist (T) | [sin(qR_{Y}), 0, cos(qR_{Y})] | K_{22}/2 |

Bend (B) | $[qRZ,0,1\u2212qRX]/q2RZ2+(1\u2212qRX)2$ | K_{33}/2 |

Δ deformation (Δ) | $[qRX,\u2212qRY,1]/q2(RX2+RY2)+1$ | 2K_{24} |

Double twist (DT) | $[\u2212qRY,qRX,1]/q2(RX2+RY2)+1$ | 2(K_{22} − K_{24}) |

Generalized | $[qRX,qRY,1+qwRZ]/q2(RX2+RY2)+(1+qwRZ)2$ | 2[K_{11} + (2 − w)K_{13} − K_{24}]a |

double splay (DS) |

^{a}

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

### B. Microscopic free energy density of deformation

^{43}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

**and with orientation $u\u0302$. For the general case of a nonuniform director field $n\u0302d$, we will use the following ansatz for the form of the ODF:**

*R**c*is a parameter that quantifies the orienting strength and $Q=\u222bdu\u0302exp{c[u\u0302\u22c5n\u0302d(R,q)]2}$ is the orientational partition function. Equation (5) 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 $[u\u0302\u22c5n\u0302d]$ 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.

*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, $adid$, is obtained by substituting the general form of the ODF [Eq. (5)] in Eq. (7); however, since the latter involves the integration over all possible orientations and the ODF is evaluated at a single point in space, $adid=auid$ is obtained.

^{45}Parsons–Lee

^{46,47}correction. In particular, the excess free energy density in

**=**

*R***0**for the deformations defined in Table I reads

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

^{48}which in turn is defined as

*e*

_{AB}= exp{−

*U*

_{AB}/

*k*

_{B}

*T*} − 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

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

^{45}In the undeformed state, the ODFs

*f*

_{d}in Eq. (8) are replaced by the corresponding position independent ODFs,

*f*

_{u}, defined in Eq. (6).

*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\u27e8uZ2\u27e9u\u22121)/2$, where ⟨⋯⟩

_{u}indicates the orientational average over the undeformed ODF defined in Eq. (6). From the profile of $\Delta addef(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,

^{49}whereas from the curvature of $\Delta addef(0,q)$ for

*q*→ 0, we can determine the elastic constants (see the expressions of $Cd$ in Table I).

### C. Explicit expressions for the elastic constants

*q*= 0 of the ODFs appearing in Eq. (8). If we introduce the following shorthand notations:

*∂*

_{q}indicates differentiation with respect to

*q*. If we introduce Eq. (11) into the expression of the microscopic free energy density [Eq. (10)], we obtain the following correspondence between the combinations of elastic constants $Cd$ and the second-order term of the expansion:

*q*= 0 has the general form

*q*= 0 has the general form

*q*= 0 are collected in Table II. We note that, since all terms involve at least one Cartesian component of

**and particle**

*R**A*is placed at the origin,

*R*_{A}=

**0**, only the second term between braces in Eq. (12) survives, so we can write

*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:

*K*

_{22}and

*K*

_{33},

*K*

_{24}:

*w*-dependent quantities in the generalized double splay deformation, we obtain the following expression for

*K*

_{13}:

^{6}relationship between

*K*

_{13}and other elastic constants is recovered,

d . | $\u2202q[gd(u\u0302,R,q)]q=02c$ . | $\u2202q\u2202q[gd(u\u0302,R,q)]q=02c$ . |
---|---|---|

S | R_{X}u_{X}u_{Z} | $RX2(uX2\u2212uZ2)\u22122RXRZuXuZ$ |

T | R_{Y}u_{X}u_{Z} | $RY2(uX2\u2212uZ2)$ |

B | R_{Z}u_{X}u_{Z} | $RZ2(uX2\u2212uZ2)+2RXRZuXuZ$ |

Δ | R_{X}u_{X}u_{Z} − R_{Y}u_{Y}u_{Z} | $RX2(uX2\u2212uZ2)+RY2(uY2\u2212uZ2)\u22122RXRYuXuY$ |

DT | R_{X}u_{Y}u_{Z} − R_{Y}u_{X}u_{Z} | $RX2(uY2\u2212uZ2)+RY2(uX2\u2212uZ2)\u22122RXRYuXuY$ |

DS | R_{X}u_{X}u_{Z} + R_{Y}u_{Y}u_{Z} | $RX2(uX2\u2212uZ2)+RY2(uY2\u2212uZ2)+2RXRYuXuY$ |

$\u22122w(RXRZuXuZ+RYRZuYuZ)$ |

d . | $\u2202q[gd(u\u0302,R,q)]q=02c$ . | $\u2202q\u2202q[gd(u\u0302,R,q)]q=02c$ . |
---|---|---|

S | R_{X}u_{X}u_{Z} | $RX2(uX2\u2212uZ2)\u22122RXRZuXuZ$ |

T | R_{Y}u_{X}u_{Z} | $RY2(uX2\u2212uZ2)$ |

B | R_{Z}u_{X}u_{Z} | $RZ2(uX2\u2212uZ2)+2RXRZuXuZ$ |

Δ | R_{X}u_{X}u_{Z} − R_{Y}u_{Y}u_{Z} | $RX2(uX2\u2212uZ2)+RY2(uY2\u2212uZ2)\u22122RXRYuXuY$ |

DT | R_{X}u_{Y}u_{Z} − R_{Y}u_{X}u_{Z} | $RX2(uY2\u2212uZ2)+RY2(uX2\u2212uZ2)\u22122RXRYuXuY$ |

DS | R_{X}u_{X}u_{Z} + R_{Y}u_{Y}u_{Z} | $RX2(uX2\u2212uZ2)+RY2(uY2\u2212uZ2)+2RXRYuXuY$ |

$\u22122w(RXRZuXuZ+RYRZuYuZ)$ |

## III. NUMERICAL IMPLEMENTATION

### A. Model details

*M*= 6, 12, and 24 spheres of diameter

*σ*. The potential energy between two rods

*A*and

*B*is defined as

*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*_{0} of a particle, used to calculate the packing fraction *ϕ* = *ρv*_{0} of the system, is simply the sum of the volumes of the hard spheres: *v*_{0} = *Mπσ*^{3}/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;^{50} a density of vertices on the surface equal to 1000 *σ*^{−2} was adopted. The results for the three particle types are reported in Table III.

### 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 $adex(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^{12}.

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, $u\u0302A$ and $u\u0302B$. 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}).

^{51}

### 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, $au=auex+auid$. 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^{52} implemented in Matlab^{53} 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.

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

*ρ*

_{iso}and the nematic

*ρ*

_{nem}phases in the coexistence regime are obtained from the solution of the following system of equations:

*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*_{1}⟩, 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*^{2}*π*. Calculations are performed for fixed values of the *c* and *q* parameters (*q* = 0 for the undeformed system).

## 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*^{2} > 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 $\Delta addef(0,q)$ for splay, twist, bend, Δ, and double twist deformations, whose curvatures provide *K*_{11} = 1.54 *k*_{B}*Tσ*^{−1}, *K*_{22} = 0.60 *k*_{B}*Tσ*^{−1}, *K*_{33} = 20.18 *k*_{B}*Tσ*^{−1}, and *K*_{24} = 2.17 *k*_{B}*Tσ*^{−1}.

Interestingly, *K*_{24} is larger than *K*_{22}, which is confirmed by direct calculation of the difference *K*_{22} − *K*_{24} 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 $CDS(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*_{13} = 3.29 *k*_{B}*Tσ*^{−1}. Thus, the saddle–splay constant turns out to be higher than all the other elastic coefficients but *K*_{33}. 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*_{24} > *K*_{22} and *K*_{24} > *K*_{11}, respectively (in the latter case there is no *K*_{13} 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*_{11} − *K*_{24}).

### 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*_{13}, 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.

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*_{33} and obey *K*_{22} < *K*_{11}, 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*_{11} and *K*_{22}, irrespective of the aspect ratio of the particles.

^{54}hence, they must be accompanied by defects or by other deformations, which give a positive contribution to the total free energy.

^{55}We must also mention that the value that we calculate for

*K*

_{11}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,

*l*is the average longitudinal projection of a particle and

*ρ*

_{0}=

*ρ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.

^{58}The dashed lines in the plots of Fig. 7 show the “renormalized” splay elastic constant of THSRs, calculated assuming $l=(M\u22121)\sigma \u27e8|uZ|\u27e9u$, 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*_{24}. 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.^{42}For the same systems, we find also*K*_{24}>*K*_{11}, 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*_{13}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*_{24}>*K*_{22}and large*K*_{13}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*_{24} > *K*_{22}, 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

*V*contained in a nematic sample of volume $V$ to the following microscopic expression for the deformation free energy:

*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. (A1) is simply the free energy of the region when the nematic director field is uniform.

*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

**and**

*R*

*R*_{A}is found, and by dropping the now superfluous

*A*subscript one can write

**.**

*R*^{32}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. (A1), one could obtain different expressions for $\Delta addef$ that in turn would lead to different values for the surface-like elastic constants. As an example, the following change of variables was proposed:

*R*_{0}-dependent

*even in absence of deformation*, through its domain of integration. Analogous conclusions were derived, from a different perspective, in Ref. 33.

^{59}

### APPENDIX B: EQUIVALENCES BETWEEN INTEGRALS

*Z*axis of the reference frame, the quantities involved in the last expression transform according to the following relations:

*Z*axis of the reference frame, the quantities involved in this expression transform according to the following relations:

*Z*axis of the reference frame, Eq. (B5) is recovered.

## REFERENCES

*K*_{24}

*K*_{13}divergence term

*K*_{24}in nematic liquid crystals

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

*An Elementary Course on the Continuum Theory for Nematic Liquid Crystals*

_{13}surfacelike elastic term in the free energy of a nematic liquid crystal

*The Physics of Liquid Crystals*

*K*

_{13}of a nematic liquid crystal

_{24}of nematic liquid crystals from an algebraic approach

*Numerical Recipes: The Art of Scientific Computing*

*Polymer Liquid Crystals*

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.