Perturbative density functional methods for cholesteric liquid crystals

TL;DR

This paper develops a numerical framework to predict the macroscopic properties of cholesteric liquid crystals from microscopic models, demonstrating high accuracy and efficiency in calculating equilibrium structures using perturbative methods.

Contribution

It introduces a perturbative approach based on Poniewierski-Stecki equations for cholesteric phases, improving computational efficiency and accuracy over previous methods.

Findings

Predictions align well with full-functional descriptions.

Perturbative treatment valid for pitch lengths of a few dozen particle diameters.

Full angle-dependent virial coefficients enhance numerical efficiency.

Abstract

We introduce a comprehensive numerical framework to generically infer the emergent macroscopic properties of uniaxial nematic and cholesteric phases from that of their microscopic constituent mesogens. This approach, based on the full numerical resolution of the Poniewierski-Stecki equations in the weak chirality limit, may expediently handle a wide range of particle models through the use of Monte-Carlo sampling for all virial-type integrals. Its predictions in terms of equilibrium cholesteric structures are found to be in excellent agreement with previous full-functional descriptions, thereby demonstrating the quantitative validity of the perturbative treatment of chirality for pitch lengths as short as a few dozen particle diameters. Furthermore, the use of the full angle-dependent virial coefficients in the Onsager-Parsons-Lee formalism increases its numerical efficiency by several orders of magnitude over that of these previous methods. The comparison of our results with numerical simulations however reveals some shortcomings of the Parsons-Lee approximation for systems of strongly non-convex particles, notwithstanding the accurate inclusion of their full effective molecular volume. Further potential limitations of our theory in terms of phase symmetry assumptions are also examined, and prospective directions for future improvements discussed.

Tables

Table 1. Table 1: Comparison of this work, Ref.  Dussi  et al. , 2015 and Ref.  Frezza  et al. , 2014 for hard helices with N s = 15 subscript 𝑁 𝑠 15 N_{s}=15 and l = 10 ​ σ 𝑙 10 𝜎 l=10\,\sigma at the isotropic-cholesteric transition. | Δ ​ 𝒫 | / 𝒫 eq Δ 𝒫 subscript 𝒫 eq \left|\Delta\mathcal{P}\right|/\mathcal{P}_{\rm eq} represents the relative difference between the equilibrium pitch predicted by our method and those reported in the references above, normalised by our value. Our results were obtained by averaging over 8 independent runs with N MC = 5 × 10 11 subscript 𝑁 MC 5 superscript 10 11 N_{\rm MC}=5\times 10^{11} ; numerical dispersion was found to be below 5% for all computed pitches, and smaller than the last significant digit for all other variables.

		$r_{.2}{p}_2$	$r_{.2}{p}_4$	$r_{.2}{p}_8$	$r_{.4}{p}_2$	$r_{.4}{p}_4$	$r_{.4}{p}_8$
$\eta$ ($\%$)	here	30.7	28.1	26.5	41.2	35.0	29.2
	Ref. Dussi et al.,2015	30.0	27.4	25.8	40.3	34.0	28.2
$S$	here	0.704	0.684	0.696	0.629	0.632	0.621
	Ref. Dussi et al.,2015	0.699	0.677	0.690	0.612	0.622	0.619
	Ref. Frezza et al.,2014	0.64	0.66	0.68	0.60	0.61	0.61
$-{10}^4\times k_t$ ($k_{b}T/{\sigma }^2$)	here	4.12	52.1	46.0	-11.6	130	139
	Ref. Dussi et al.,2015	4.27	47.3	42.3	-10.7	110	115
	Ref. Frezza et al.,2014	4.83	41.36	29.03	-3.83	98.35	110.13
$K_2$ ($k_{b}T/\sigma$)	here	0.201	0.197	0.207	0.170	0.167	0.181
	Ref. Dussi et al.,2015	0.199	0.194	0.203	0.160	0.150	0.136
	Ref. Frezza et al.,2014	0.154	0.177	0.184	0.153	0.152	0.159
${𝒫}_{\mathrm{eq}}$ ($\sigma$)	here	-3065	-238	-283	921	-81	-82
	Ref. Dussi et al.,2015	-2990	-260	-310	965	-93	-90
	Ref. Frezza et al.,2014	-2008	-268	-399	2509	-97	-90
$\left|\Delta 𝒫\right|/{𝒫}_{\mathrm{eq}}$ ($\%$)	Ref. Dussi et al.,2015	2.45	9.24	9.54	4.78	14.8	9.76
	Ref. Frezza et al.,2014	34.5	12.6	41.0	172	19.8	9.76

Equations

$$\begin{split}\mathscr{F}=&\>\mathscr{F}_{0}+\frac{1}{2}\int_{V}d\mathbf{r}\times\Big{\{}K_{1}\left(\nabla\cdot\mathbf{n}\right)^{2}+K_{2}\left(\mathbf{n}\cdot\left[\nabla\times\mathbf{n}\right]\right)^{2}\\ &+K_{3}\left(\mathbf{n}\times\left[\nabla\times\mathbf{n}\right]\right)^{2}+2k_{t}\left(\mathbf{n}\cdot\left[\nabla\times\mathbf{n}\right]\right)\Big{\}}\end{split}$$

$$\rho (\mathbf{r},\mathcal{R})=\rho (\mathbf{r},\mathbf{u}).$$

$$\oint d\mathcal{R}\times \rho (\mathbf{r},\mathbf{u})=2\pi \oint d\mathbf{u}\times \rho (\mathbf{r},\mathbf{u})=\frac{N}{V}\forall \mathbf{r}\in V$$

$$\rho (\mathbf{r},\mathbf{u})=\rho \times \psi \left[\mathbf{n}(\mathbf{r})\cdot \mathbf{u}\right]$$

$$\mathbf{n}(\mathbf{r})={\mathbf{n}}_0\equiv {\mathbf{e}}_z\forall \mathbf{r}\in V.$$

$$\rho (\mathbf{r},\mathbf{u})=\rho (\mathbf{u})=\rho \times \psi \left(\cos \theta \right)$$

$${\mathcal{F}}_0^{\infty }[\rho ]={\mathcal{F}}_0^{\infty }[\psi ]={\mathcal{F}}_{\mathrm{id}}[\psi ]+{\mathcal{F}}_{\mathrm{ex}}^{\infty }[\psi ],$$

$$\begin{split}\frac{\beta\mathscr{F}_{\rm id}[\psi]}{V}=&\>\rho(\log\rho-1)+4\pi^{2}\rho\int_{0}^{\pi}d\theta\\ &\times\sin\theta\log\Big{\{}\psi(\cos\theta)\Big{\}}\psi(\cos\theta),\end{split}$$

$$\begin{array}{rl}\frac{\beta {\mathcal{F}}_{\mathrm{ex}}^{\infty }[\psi ]}{V}=& \iint \\ & \times f({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2)\psi ({\mathbf{n}}_1\cdot {\mathbf{u}}_1)\psi ({\mathbf{n}}_2\cdot {\mathbf{u}}_2)\end{array}$$

$${\mathbf{n}}_i\cdot {\mathbf{u}}_i\equiv \mathbf{n}({\mathbf{r}}_i)\cdot {\mathbf{u}}_i=\cos {\theta }_i\forall i\in [[1,2]].$$

$$f(\mathbf{r}_{12},\mathcal{R}_{1},\mathcal{R}_{2})=\exp\Big{\{}-\beta U(\mathbf{r}_{12},\mathcal{R}_{1},\mathcal{R}_{2})\Big{\}}-1$$

$$\begin{split}\psi^{(\rho)}_{\rm eq}(\cos\theta)=&\>\frac{1}{Z}\times\exp\bigg{\{}-\frac{\rho}{4\pi^{2}}\int_{0}^{\pi}d\theta^{\prime}\\ &\times\sin\theta^{\prime}\frac{E^{\infty}(\theta,\theta^{\prime})+E^{\infty}(\theta^{\prime},\theta)}{2}\psi^{(\rho)}_{\rm eq}(\cos\theta^{\prime})\bigg{\}}\end{split}$$

$$\begin{array}{rl}{E}^{\infty }({\theta }_1,{\theta }_2)=& -{\int }_{V}d{\mathbf{r}}_{12}{\iint }_0^{2\pi }d{\alpha }_{1}d{\alpha }_2{\iint }_0^{2\pi }d{\varphi }_{1}d{\varphi }_2\\ & \times f({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2)\end{array}$$

$$Z=4{\pi }^2{\int }_0^{\pi }d\theta \times \sin \theta \times {\psi }_{\mathrm{eq}}^{(\rho )}(\cos \theta ).$$

$$\frac{\delta^{2}\beta\mathscr{F}_{\rm ex}[\psi]}{\delta\psi(\mathbf{n}_{1}\cdot\mathbf{u}_{1})\delta\psi(\mathbf{n}_{2}\cdot\mathbf{u}_{2})}\Biggr{|}_{\rho\psi}=-\rho^{2}\times c^{(2)}(\mathbf{r}_{12},\mathcal{R}_{1},\mathcal{R}_{2};[\rho\psi])$$

$$\underset{\rho \to 0}{\lim }{c}^{(2)}({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2;[\rho \psi ])=f({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2).$$

$$\begin{array}{rl}\frac{\beta {\mathcal{F}}_{\mathrm{ex}}[\psi ]}{V}=& \iint \\ & \times \overline{c^{(2)}}({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2;[\rho \psi ])\psi ({\mathbf{n}}_1\cdot {\mathbf{u}}_1)\psi ({\mathbf{n}}_2\cdot {\mathbf{u}}_2)\end{array}$$

$$\begin{array}{rl}\overline{c^{(2)}}({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2;[\rho \psi ])=& 2{\int }_0^{1}d\lambda {\int }_0^{\lambda }d{\lambda }^{\prime }\\ & \times c^{(2)}({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2;[{\lambda }^{\prime }\times \rho \psi ]).\end{array}$$

$$\overline{c^{(2)}}({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2;[\rho \psi ])\simeq f({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2)$$

$$c^{(2)}({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2;[\rho \psi ])=c_{\mathrm{ref}}^{(2)}\left\{\frac{\Vert {\mathbf{r}}_{12}\Vert }{\sigma ({\mathbf{u}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2)};\rho \right\}$$

$$\overline{c^{(2)}}({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2;[\rho \psi ])=G(\rho v_{\mathrm{ref}})\times f({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2)$$

$${\eta }_{\mathrm{ref}}^{(\rho )}\equiv \rho v_{\mathrm{ref}}={\eta }^{(\rho )}\times \frac{v_{\mathrm{ref}}}{v}.$$

$$G_{\mathrm{PW}}(\eta )$$

$$G_{\mathrm{PL}}(\eta )$$

$$\underset{\eta \to 0}{\lim }G(\eta )=1.$$

$${\mathcal{F}}_{\mathrm{ex}}[\psi ]$$

$$E({\theta }_1,{\theta }_2)$$

$${\mathcal{F}}_0[\psi ]={\mathcal{F}}_{\mathrm{id}}[\psi ]+G_{\mathrm{PL}}\left({\eta }_{\mathrm{ref}}^{(\rho )}\right)\times {\mathcal{F}}_{\mathrm{ex}}^{\infty }[\psi ]$$

$$\begin{split}\psi^{(\rho)}_{\rm eq}(\cos\theta)=&\>\frac{1}{Z}\times\exp\bigg{\{}-G_{\rm PL}\left(\eta_{\rm ref}^{(\rho)}\right)\times\frac{\rho}{4\pi^{2}}\int_{0}^{\pi}d\theta^{\prime}\\ &\times\sin\theta^{\prime}\frac{E^{\infty}(\theta,\theta^{\prime})+E^{\infty}(\theta^{\prime},\theta)}{2}\psi^{(\rho)}_{\rm eq}(\cos\theta^{\prime})\bigg{\}}.\end{split}$$

$$\overline{c^{(2)}}({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2;[\rho \psi ])\simeq c^{(2)}({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2;[\rho \psi ]).$$

$$\begin{array}{rl}\beta K_1[\psi ]=& \iint \\ & \times c^{(2)}({\mathbf{r}}_{12},{\mathcal{R}}_1,{\mathcal{R}}_2;[\rho \psi ])\dot{\psi }(u_{1z})\dot{\psi }(u_{2z})r_x^2{u}_{1x}{u}_{2x}\end{array}$$

Full text

I Introduction

Liquid crystals (LCs) are a fascinating example of self-assembling soft-matter systems, characterised by the singular coexistence of long-range order and fluidity. These organised phases typically arise from a specific set of structural and chemical features displayed by their constituent particles,de Gennes and Prost (1993) the most notable of which being shape anisotropy. In other words, these particles are elongated along one or more directions, and therefore generally possess both a long and a short axis — respectively defined as the directions of maximum and minimum extent of their molecular backbone.

A remarkable consequence of this microscopic anisotropy is the wealth of symmetry-breaking transitions involved in LC phase-ordering kinetics, as not only the positions, but also the orientations of the molecules may display long-range correlations.de Gennes and Prost (1993) Therefore, LCs are characterised by a degree of orientational order, potentially associated with a form of positional organisation — leading to an impressive variety of macroscopic structures.Goodby and Gray (2008) Understanding the link between this range of phase symmetries and the underlying properties of their constituent particles has remained a long-standing challenge of soft-condensed matter theory.Wilson (2005)

The simplest example of LC assembly is the so-called uniaxial nematic phase, in which one of the molecular axes points on average in a particular direction called the nematic director, denoted by $\mathbf{n}$, while retaining fluid-like positional disorder. The resulting phases may then be qualified as prolate or oblate, respectively referring to the local orientational order of the particle long or short axes. This structural anisotropy bestows upon nematic LCs a strong optical birefringence,Wu (1986) much like regular uniaxial crystals, along with unique electricalFréedericksz and Zolina (1933) and rheologicalDoi (1981) properties. In the absence of external aligning forces, the director generally fluctuates throughout the sample, and is therefore usually referred to as a director field $\mathbf{n}(\mathbf{r})$, with $\mathbf{r}$ the spatial coordinate vector.

An intriguing special case of this spatial modulation corresponds to the cholesteric phase, typically formed by systems of chiral particles, in which $\mathbf{n}(\mathbf{r})$ periodically rotates along an axis normal to the local director, as illustrated in Fig. 1. Hence, cholesteric LCs are commonly described in terms of a single lengthscale $\mathcal{P}$ — termed the cholesteric pitch — defined as the spatial period of the director rotation, and fully characterising their macroscopic structure.

As the earliest LC phase to be identified experimentally in 1888,Reinitzer (1888) cholesterics boast a remarkably long history along with a rich array of applications; the most notable example indubitably lies in the twisted nematic cell, which single-handedly triggered the advent of modern LC display technology.Dunmur and Sluckin (2011) However, an impressive variety of further uses of cholesteric systems have been proposed in recent years, ranging from physicalMoreira et al. (2004) and chemicalMujahid et al. (2010) sensors to surface-morphing applications,Liu and Broer (2015) including the design of tunable-wavelength lasersFinkelmann et al. (2001) and core-sheath fibers for smart textiles.Enz and Lagerwall (2010)

Their original observation in cholesterol derivatives marked the discovery of the thermotropic class of LCs, in which phase transitions are governed by temperature variations. However, lyotropic cholesteric phases, in which transitions are mostly ruled by mesogen concentration, are also commonly observed in many bio-colloidalRobinson (1961); Livolant and Leforestier (1996); Dogic and Fraden (2006); Uematsu and Uematsu (1984); Hiltrop (2001); Sato and Teramoto (1996); Werbowyj and Gray (1976); Revol and Marchessault (1993) and polymerAharoni (1979); Sato et al. (1993); Watanabe, Kamee, and Fujiki (2001) solutions. A fundamental feature of both systems is the high sensitivity of the pitch $\mathcal{P}$ to the chemical structure and thermodynamic state of their constituent particles, which provides the basis for many of the above applications; the nature of this intricate dependence has been intensively investigated experimentally over the last decades in a number of different contexts.Yevdokimov, Skuridin, and Salyanov (1988); Van Winkle et al. (1990); Stanley, Hong, and Strey (2005); Dogic and Fraden (2000); Grelet and Fraden (2003); DuPré and Duke (1975); Dong and Gray (1997); Miller and Donald (2003)

Conversely, theoretical and numerical descriptions of the cholesteric phase have long been hampered by both the complexity of the microscopic chiral interactions and the multi-scale nature of the problem, as cholesteric pitches measured experimentally are usually several orders of magnitudes larger than typical particle sizes.de Gennes and Prost (1993) This large spatial periodicity renders direct atomistic simulations largely impractical, due to both the substantial size of the simulation boxes required and the unsuitability of standard periodic boundary conditions.Frenkel (2013) Furthermore, the quick escalation of computational costs with increasing particle aspect ratios largely restricts their applicability to systems of short, highly chiral particles — in stark contrast with the properties of many experimentally-relevant cholesteric mesogens.Van Winkle et al. (1990); Livolant and Leforestier (1996); Stanley, Hong, and Strey (2005); Dogic and Fraden (2006, 2000); Grelet and Fraden (2003); Usov et al. (2015); Wang, Hamad, and MacLachlan (2016)

The majority of theoretical attempts to overcome these limitations have thus chosen to follow in the footsteps of Onsager,Onsager (1949) whose seminal microscopic description of the nematic state provides a natural framework to investigate the self-assembly of highly anisotropic particles. This original approach, which prefigures the modern formulation of classical density functional theory (DFT),Evans (1992a) has been successfully used to predict the isotropic-nematic (I-N) phase transition of a number of experimental systems,Fraden et al. (1989); Lagerwall and Scalia (2008); Buining and Lekkerkerker (1993); Zhang and van Duijneveldt (2006) and was subsequently generalised to the cholesteric phase by StraleyStraley (1976) in the limit of weak chirality.

However, most implementations of Straley’s theory have so far been limited by stringent analytical assumptions to highly simplified particle models, in which a number of further ad hoc approximations had to be introduced for the tractable inclusion of chiral contributions in the virial formalism of the extended Onsager theory.Osipov (1985); Odijk (1987); Pelcovits (1996); Evans (1992b); Varga and Jackson (2006, 2011); Wensink and Jackson (2009, 2011); Wensink (2014) Some numerical attempts to tackle more realistic systems — including both steric and soft repulsion — represented significant steps forward, but were still hindered by a cumbersome mathematical formalism and numerical procedure.Tombolato and Ferrarini (2005); Tombolato, Ferrarini, and Grelet (2006); Frezza et al. (2014) In particular, their recourse to the partial decomposition of the particle density field in a basis of rotational invariants significantly limited both their computational efficiency and numerical precision due to potential finite-order truncation effects.Lekkerkerker et al. (1984)

A notable exception to the previous pattern lies in the recent works of Belli et al.,Belli et al. (2014); Dussi et al. (2015) in which a unified framework capable of handling arbitrarily complex particle shapes and interaction potentials independent of Straley’s perturbative description was introduced. This approach, based on the use of Monte-Carlo (MC) integration for the computation of all generalised virial coefficients, however suffers from two major drawbacks. Firstly, it requires the evaluation of the full chiral free energy landscape of the molecules using a finite grid of macroscopic pitches, leading to high computational costs along with discretisation effects stemming from the limited grid resolution. Secondly, the reliable calculation of equilibrium pitches typically requires several preliminary runs in order to narrow down the pitch grid to the relevant parameter range, further increasing its computational expense.

On these premises, we here aim to make progress towards an efficient, quantitatively reliable theoretical description of the cholesteric structure and underlying nematic properties of systems of arbitrary particles. Combining the numerical integration techniques of Refs. Belli et al., 2014; Dussi et al., 2015 with Straley’s perturbative frameworkStraley (1976) enables us to circumvent the limitations of previous approaches, and allows for the expedient computation of all relevant variables with an adjustable level of precision for a wide variety of particle models. Furthermore, the relative simplicity of the corresponding mathematical expressions lends itself to a thorough examination of the different physical assumptions underpinning the DFT formalism, and of their potential shortcomings in the case of chiral nematic phases.

The structure of the paper is thus organised as follows. In order that the origin of the potential shortcomings of the theory can be fully understood, we first provide in Sec. II a comprehensive account of the theory and numerical implementation, focussing on the physical nature of the approximations that are being made. Specifically, we outline in Sec. II.1 the original Onsager DFT for nematic LCs, with a particular emphasis on phase symmetry assumptions, and focus in Sec. II.2 on the approximations relative to the treatment of correlations between finite-size particles. We tackle in Sec. II.3 the applications of the theory to the microscopic determination of LC elasticities, which are intimately linked to the perturbative treatment of chirality in cholesterics, before detailing the computational procedure used in Sec. II.4. We then devote Sec. III to the extensive comparison of our results to previous numerical and theoretical studies of several model systems, and discuss potential deficiencies of the different density-functional assumptions stemming from both particle non-convexity and phase biaxiality. We finally conclude in Sec. IV with a few remarks on the possible consequences of these shortcomings on the general reliability of our approach, and highlight some prospective directions for future research.

II Microscopic theory of liquid-crystal elasticities

We here choose to model cholesterics as weakly distorted uniaxial nematic phases, whose bulk free energies are typically described at the continuum level by the standard Oseen-Frank field functionalFrank (1958)

[TABLE]

with $k_{b}T$ the thermal energy and $V$ the total volume of the sample. $\mathscr{F}_{0}$ corresponds to the free energy density of the uniformly aligned nematic phase, and $k_{t}$ is a chiral strength parameter which vanishes in the case of achiral mesogens. $K_{1}$, $K_{2}$ and $K_{3}$ are referred to as the Frank elastic constants, and quantify the respective free energy penalties associated with splay, twist and bend deformation modes of the local director field.de Gennes and Prost (1993) It should be noted that this compact expression derives from a finite-order truncation of the free energy gradient expansion, and is therefore only valid in the limit of long-wavelength director fluctuations — an assumption consistent with the characteristically large values of experimental cholesteric pitches.Yevdokimov, Skuridin, and Salyanov (1988); Van Winkle et al. (1990); Stanley, Hong, and Strey (2005); Dogic and Fraden (2000); Grelet and Fraden (2003); DuPré and Duke (1975); Dong and Gray (1997); Miller and Donald (2003)

Eq. (1) thus provides a natural theoretical framework for the description of smoothly modulated nematic phases, and has been successfully used in recent years for the analyticalMajumdar and Lewis (2017); Lewis et al. (2014) and numericalMajumdar and Lewis (2017); Dammone et al. (2012) modelling of a number of experimentally-relevant systems. However, its main shortcoming lies in the coarse-graining of all bulk particle properties into the three $K_{i}$ scalars — understood as free parameters arising from basic mean-field symmetries, and often taken to be equal for the sake of simplicity.de Gennes and Prost (1993) In order to provide a full quantitative, bottom-up description of nematic liquid crystals, it is therefore desirable to derive microscopic expressions for these Frank elastic constants, and thus provide a link between mesogen characteristics and phase properties in such self-assembled systems.

A typical starting point for this analysis, as first proposed in Straley’s pioneering work,Straley (1976) lies in the molecular description of the free energy $\mathscr{F}_{0}$ of the reference uniform nematic state. The first tractable expression for this free energy was reported by OnsagerOnsager (1949) for generic systems of highly anisotropic particles and forms the basis of his classical density functional theory, which we outline in the next section.

II.1 The Onsager DFT for uniaxial LCs

The original formulation of the Onsager DFT aims to quantify the onset of macroscopic nematic order in terms of competing entropic contributions at the microscopic level.Onsager (1949) Similarly to other DFT approaches,Evans et al. (2016) its starting point lies in the derivation of the free energy of the system as a functional of the mesoscopically-averaged single-particle density $\rho(\mathbf{r},\mathcal{R})$, defined as the number of particles at position $\mathbf{r}$ with orientation $\mathcal{R}$ per unit volume and solid angle. $\mathcal{R}$ is generally represented by a full $3\times 3$ rotation matrix, which we choose to parametrise in the $z$-$y$-$z$ convention using the three Euler angles $\alpha$, $\theta$ and $\phi$, as depicted in Fig. 1a. In the special case of uniaxial calamitic particles, the internal rotational degree of freedom $\alpha$ becomes irrelevant, and the symmetries of the system are fully captured by the respective polar and azimuthal angles $\theta$ and $\phi$, defining the orientation vector $\mathbf{u}$ of the particle long axis.

Analogously, in the case of prolate uniaxial phases, internal particle angles $\alpha$ do not display long-range correlations, and the macroscopic structure of the sample is thus entirely characterised by the spatial distribution of the long axes of its constituent particles. The corresponding density field hence reduces to

[TABLE]

Furthermore, owing to the lack of positional order in nematic LCs, the number of particles per unit volume is often taken to be spatially uniform within the sample, which yields the following marginal normalisation constraint:

[TABLE]

with $N$ the total number of particles in the system. Denoting this uniform number density by $\rho\equiv N/V$, the density field of a bulk nematic phase with director $\mathbf{n}$ then takes the general formWensink and Jackson (2009)

[TABLE]

with $\psi\left[\mathbf{n}(\mathbf{r})\cdot\mathbf{u}\right]$ the so-called orientation distribution function (ODF), which quantifies the probability of finding a particle at point $\mathbf{r}$ with its long axis borne by $\mathbf{u}$ with respect to the local director $\mathbf{n}(\mathbf{r})$. Note that ansatz (4) further reduces the orientational dependence of the density field to the single scalar $\mathbf{n}(\mathbf{r})\cdot\mathbf{u}$, and thereby assumes that the system retains local cylindrical invariance.Straley (1976) This symmetry constraint thus restricts its applicability to the description of weakly distorted nematic phases, consistent with the Frank-Oseen perturbative framework introduced above.

Following in the footsteps of Onsager, let us first consider the case of a spatially uniform director field $\mathbf{n}_{0}$, which we define as the $z$-axis of our reference frame

[TABLE]

Eq. (4) then reduces to

[TABLE]

with $\theta$ the particle polar angle as defined in Fig. 1a. The free energy of the system can then be cast in the usual formOnsager (1949)

[TABLE]

where the subscript [math] denotes the total free energy of the undeformed phase, and $\mathscr{F}_{\rm id}$ and $\mathscr{F}_{\rm ex}^{\infty}$ its respective ideal and excess components,

[TABLE]

[TABLE]

with $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ the respective positions of the particle centers of mass, $\mathbf{r}_{12}\equiv\mathbf{r}_{2}-\mathbf{r}_{1}$ and

[TABLE]

The ideal free energy (Eq. (8)) corresponds to the free energy of an ideal gas of anisotropic particles with uniform density $\rho$, whose two terms can be interpreted as the respective contributions of translational and rotational entropy. The excess free energy (Eq. (9)) stems from the truncation of the Mayer cluster expansionMayer and Mayer (1940) at the second-virial level, which as discussed in the next section is only exact in the limit of infinitely thin particles,Onsager (1949) and is hence denoted by an $\infty$ superscript. It accounts for microscopic interactions in the pairwise additive approximation through the so-called Mayer $f$-functionMayer and Mayer (1940)

[TABLE]

with $\beta$ the inverse thermal energy and $U$ the particle pairwise potential.

It can then be proven that the equilibrium ODF $\psi^{(\rho)}_{\rm eq}$ that minimises the nematic free energy (Eq. (7)) at fixed density $\rho$ satisfies the following self-consistent equation:

[TABLE]

with $E$ the orientation-dependent second-virial coefficientBelli et al. (2014)

[TABLE]

and $Z$ a normalisation constant ensuring that Eq. (3) is satisfied, namely,

[TABLE]

Note that in the case of monodisperse systems of identical molecules, it follows from the invariance of the interaction potential with respect to particle exchange that $E^{\infty}(\theta_{1},\theta_{2})=E^{\infty}(\theta_{2},\theta_{1})$. We however choose to keep the explicit symmetric dependence of Eq. (12) on $E^{\infty}$ for reasons of numerical convergence, as this form provides an effective average over the non-diagonal coefficients of $E^{\infty}$ obtained using the computational procedure introduced in Sec. II.4.

Eqs. (6) and (12) thus fully describe the equilibrium structure of a nematic LC at any given density $\rho$. However, this original form of the Onsager DFT is restricted by geometrical assumptions to infinitely thin particles and uniform director fields, and is therefore of limited practical relevance. We hence devote the next two sections to a thorough, though non-exhaustive, overview of subsequent attempts to overcome both these limitations.

II.2 Introducing finite mesogen anisotropy

The main difficulty in dealing with finite-size particles lies in the need to go beyond the second-virial approximation, as the relative weight of higher order contributions rapidly increases with decreasing aspect ratios;Frenkel (1987) the Onsager theory has thus been found to be quantitatively inaccurate for systems with aspect ratios of up to 100.Straley (1973a) Therefore, several improvements have been proposed in the last decades to extend this framework to a wider class of experimentally-realistic mesogens.

While some studies chose to explicitly include the third and higher order virial terms in the excess free energy expansion (Eq. (9)), the quickly escalating computational complexity of such approaches with increasing truncation order largely restricts their applicability to highly simplified particle models.Masters (2008) Furthermore, evidence of the decreased radius of convergence of the virial series for highly anisotropic mesogens may paradoxically limit their numerical stability range to regions of lower density,You, Vlasov, and Masters (2005) and even prevent the prediction of a stable nematic phase altogether for aspect ratios as low as 10 as one pushes the expansion beyond the third order in $\rho$.You et al. (2012)

A somewhat orthogonal approach is to forego this virial formalism as a starting point, and instead exploit the more general relation between excess free energy and inter-particle correlations in the density functional framework,Hansen and McDonald (1976) which in our case takes the form

[TABLE]

using the notation of Sec. II.1. $c^{(2)}$ is known as the two-particle direct correlation function (DCF), and corresponds to the “direct” component of the usual pair correlation function (PCF), omitting indirect contributions stemming from mutual correlations with neighbouring particles. In the low-density limit, this DCF can be shown to reduce to the Mayer $f$-function (Eq. (11))Vroege and Lekkerkerker (1992)

[TABLE]

However, this equality once again falters in the case of particles with a lower aspect ratio, for which the onset of liquid-crystalline order occurs at higher mesogen concentrations. In the general case, the exact relation between DCF, director profile and excess free energy at given uniform density $\rho$ can be obtained from Eq. (15) by double functional integration,

[TABLE]

with $\overline{c^{(2)}}$ the DCF integrated along a compression path from a reference state with particle density $\rho=0$,Vroege and Lekkerkerker (1992)

[TABLE]

Note that plugging Eq. (16) into Eq. (18) yields the density-functional equivalent of the second-virial approximation

[TABLE]

which reduces Eq. (17) to the Onsager excess free energy (Eq. (9)). Hence, Eqs. (17) and (18) render the accuracy of the excess free energy functional purely contingent on that of the underlying DCF for systems of arbitrarily anisotropic particles.

While good approximations for the DCF of spherical particle fluids have been obtained by directly solving the Ornstein-Zernike equation using various closure relations,Wertheim (1963); Morita and Hiroike (1960) the extra orientational degrees of freedom render this process mostly intractable when particle anisotropy is introduced.Vroege and Lekkerkerker (1992) A number of recent efforts to derive closed-form expressions for the excess free energy and associated DCF of anisotropic systems have thus resorted to first-principle geometric considerations instead, on the basis of Rosenfeld’s original fundamental measure theory (FMT) for hard spheres.Rosenfeld (1988) Although FMT presumably provides a good description of the DCF spatial-orientational coupling, such approaches are unfortunately limited by construction to convex particle modelsWittmann, Marechal, and Mecke (2016) and are not naturally suited for the treatment of soft interaction potentials.Schmidt (1999)

Another class of approaches revolve around the derivation of ad hoc approximate expressions for these quantities by remapping either the relevant DCF or excess free energy to those of the well-characterised hard-sphere fluid.Vroege and Lekkerkerker (1992) A significant number of such methods rely on the so-called decoupling approximation, in which inter-particle correlations are accounted for via simplified ansätze of the form

[TABLE]

irrespective of $\psi$, with $\lVert\cdot\rVert$ the Euclidean norm and $\mathbf{u}_{12}\equiv\mathbf{r}_{12}/\lVert\mathbf{r}_{12}\rVert$. $c^{(2)}_{\rm ref}$ then corresponds to a heuristic DCF inferred from that of a reference hard-sphere system at the same density $\rho$, either explicitly through the use of a closure relation as first proposed by PynnPynn (1974) and Wulf,Wulf (1977) or implicitly from the corresponding virial expansion as suggested independently by ParsonsParsons (1979) and Lee.Lee (1987) In this framework, particle anisotropy is introduced solely through the orientation-dependent distance of nearest approach between two particles $\sigma(\mathbf{u}_{12},\mathcal{R}_{1},\mathcal{R}_{2})$. In the case of purely hard interactions, Eqs. (17) and (18) then lead to effective expressions for $\overline{c^{(2)}}$ of the general form

[TABLE]

with $v_{\rm ref}$ the molecular volume of the reference hard-sphere fluid, which may a priori differ from that of the anisotropic particles considered, as discussed in Sec. II.4. In the general case, the effective volume fraction $\eta_{\rm ref}$ of the reference system is thus linked to that of the physical system with molecular volume $v$ at identical number density $\rho$ through the relation

[TABLE]

The density-dependent prefactor $G$ is then related to the compressibility factor of the chosen reference framework — namely, the Percus-Yevick free energyPercus and Yevick (1958); Wertheim (1963) in the Pynn-Wulf (PW) approximation, and the Carnahan-Starling equation of stateCarnahan and Starling (1969) in the Parsons-Lee (PL) approach. The respective corresponding expressions for $G$ read asVroege and Lekkerkerker (1992)

[TABLE]

which both trivially satisfy

[TABLE]

Hence, Eq. (21) simply reduces to the second-virial approximation Eq. (19) in the low-density limit.

In the following, we choose to account for finite particle anisostropy through the Parsons-Lee formalism (Eq. (24)), due to both its simplicity and well-documented success in describing the nematic behaviour of a wide range of hardCamp et al. (1996); Cuetos et al. (2007) and softCuetos et al. (2005) particle models; however, most of the methods described in the next sections are easily generalisable to other theoretical frameworks. Plugging Eq. (21) into Eq. (17) then yields the relevant excess free energy and angle-dependent virial coefficient as simple rescaled forms of their respective Onsager expressions Eqs. (9) and (13),

[TABLE]

and the equilibrium properties of the system may be fully determined by substituting Eq. (26) for $\mathscr{F}_{\rm ex}^{\infty}$ in the total free energy (Eq. (7)),

[TABLE]

and Eq. (27) for $E^{\infty}$ in the self-consistent Eq. (12),

[TABLE]

We conclude this section by remarking that computational methods for inferring accurate DCFs from the PCFs measured in direct molecular simulations have also been proposed,Phuong, Germano, and Schmid (2001) and that Eq. (17) thus conceptually enables one to calculate the “exact” excess free energy of bulk nematic phases for a wide range of particle models by simply plugging the corresponding numerical DCFs into Eq. (18). Although the appeal of such a prospect is undeniable, its practical applicability is greatly hindered by the functional dependence of the DCF on the nematic density field (Eq. (4)), as $c^{(2)}$ needs to be numerically evaluated for all density distributions along the integration path (Eq. (18)). A useful approximation may therefore be borrowed from the DFT of simple fluids by invoking linear response theory in the limit of weakly non-uniform systems, which using our notations yieldsEbner, Saam, and Stroud (1976); Saam and Ebner (1977)

[TABLE]

While the validity of the assumptions underlying Eq. (30), which in our case amount to treating the nematic phase as a weakly disturbed isotropic liquid, is clearly questionable, this relation provides a simple and direct way to check the effective correlations derived from the previous theories Eqs. (21)–(24) against the numerical DCFs predicted by simulationsPhuong, Germano, and Schmid (2002) in the vicinity of the I-N transition, and is worth keeping in mind for the assessment of future theoretical developments.

II.3 Introducing weak nematic fluctuations

The inclusion of weak director fluctuations in the previous framework as first proposed by StraleyStraley (1976) relies on the two following related assumptions, both typical of a perturbative treatment: (i) that the degree of local order is identical in the distorted and undistorted phases, so that $c^{(2)}$ and $\psi$ may be taken as equal in both systems, and (ii) that the wavevectors $\mathbf{q}$ of the director spatial fluctuations are small, so that the Frank free energy expansion Eq. (1) is well defined. On the basis of these two hypotheses, a general microscopic description of the Frank elastic constants can be straightforwardly derived within the density functional framework of the previous section; the resulting set of expressions are collectively known as the Poniewierski-Stecki equations:Poniewierski and Stecki (1979)

[TABLE]

[TABLE]

[TABLE]

which in the case of chiral particles can easily be supplemented by a fourth equation for the chiral strength parameter:

[TABLE]

with $r_{j}=\mathbf{r}_{12}\cdot\mathbf{e}_{j}$ and $u_{ij}=\mathbf{u}_{i}\cdot\mathbf{e}_{j}$. It is worth emphasising that Eqs. (31)–(34) can be single-handedly inferred from Eq. (15) using solely the two previous assumptions, and are therefore not contingent on a particular expression for the nematic free energy.

In our case, this formulation is unfortunately of little use as the Parsons-Lee theory does not provide an explicit expression for the two-particle DCF $c^{(2)}$. Therefore, a more practical derivation may be obtained by following Ref. Wensink and Jackson, 2009 and explicitly performing the Taylor expansion of the distorted free energy $\mathscr{F}$ around that of the undistorted phase (Eq. (28)). Using standard expressions for the different nematic deformation modesStraley (1973b) along with the two above assumptions enables one to readily recover the Poniewierski-Stecki equations from Eqs. (9) and (28), with the effective DCF

[TABLE]

which interestingly corresponds to the first-order term of the DCF obtained through different arguments by Somoza and Tarazona for hard spherocylinders,Somoza and Tarazona (1989a) as discussed further in Sec. III.2.

We here briefly outline this derivation in the case of a cholesteric distortion, and refer the interested reader to Ref. Wensink and Jackson, 2009 for the more detailed calculations. The director field then exhibits a helical modulation described by the simple expression

[TABLE]

where we used the conventions of Fig. 1b along with assumption (ii), ensuring that $q>0$ corresponds to a right-handed phase. Eq. (36) then trivially leads to

[TABLE]

Following the previous procedure, plugging Eq. (36) into the Onsager excess free energy (Eq. (9)) and expanding the ODF up to order 2 in $q$ allows one to recast the equilibrium free energy (Eq. (28)) in the formWensink and Jackson (2009)

[TABLE]

which by term-to-term comparison with Eq. (1) yields Eqs. (32) and (34) along with the effective DCF (Eq. (35)). Similar calculations using the corresponding distorted director fields lead to the other two elastic constants,Straley (1973b) with the notable difference that the linear terms of the associated free energy expansions vanish for non-polar mesogens.de Gennes and Prost (1993)

Assumption (i) further ensures that the equilibrium ODF $\psi^{(\rho)}_{\rm eq}$ of the cholesteric phase is still given by Eq. (29) irrespective of $q$, so that the equilibrium pitch $\mathcal{P}_{\rm eq}$ minimising the distorted free energy (Eq. (40)) at given density $\rho$ is fully determined by

[TABLE]

which quantifies the competition between the chiral torque $k_{t}$, favoring a locally-twisted phase, and the Frank twist elastic constant $K_{2}$, aiming to restore local alignment.Straley (1976)

It should be noted that this perturbative treatment of chirality — along with most discussions of this paper — is based on a uniaxial ansatz (Eq. (6)) for the nematic density field, which assumes a uniform distribution of particle angles $\alpha$ and $\phi$ within any mesoscopic sample of the system. This hypothesis thus precludes any long-range orientational order for the particle short axes, and therefore dictates that phase chirality arise solely from short-range correlations prescribed by the DCF (Eq. (35)). While consistent with our weak deformation framework, this conjecture is hard to justify a priori for systems of chiral particles and may adversely affect the reliability of the theory for cholesteric phases, as discussed in Sec. III.3. We however leave these considerations aside for now in order to facilitate the comparison of our results with previous theoretical studies, in which similar symmetry assumptions were used.

II.4 Numerical procedure

We now move on to more practical considerations, and outline the general numerical methods used to infer the equilibrium properties of chiral nematic phases from that of their arbitrary constituent molecules in our density functional framework. Following the considerations of Secs. II.2 and II.3, we hence consider from here onwards systems described by the uniaxial Onsager-Parsons-Lee equilibrium ODF $\psi_{\rm eq}^{(\rho)}$ (Eq. (29)), and introduce the shorthand

[TABLE]

for all functionals $\mathscr{X}$ of the nematic density distribution.

The working principle of our algorithm can then be summarised as follows.

1. (i)

We start by numerically integrating the angle-dependent virial coefficient $E(\theta_{1},\theta_{2})$ (Eq. (27)) for the chosen particle model using a discrete grid for polar angles $\theta\in[0,\pi[$. 2. (ii)

$E$ is subsequently plugged into Eq. (29), which is solved self-consistently using a standard iterative schemevan Roij (2005) to compute the equilibrium ODF $\psi^{(\rho)}_{\rm eq}$ for a given value of $\rho$. 3. (iii)

$\psi^{(\rho)}_{\rm eq}$ is then plugged into Eqs. (31)–(34) to work out $K_{i}(\rho)$ and $k_{t}(\rho)$ for $i\in\llbracket 1,3\rrbracket$ by further numerical integration. 4. (iv)

The equilibrium pitch $\mathcal{P}_{\rm eq}$ of the system at the chosen density $\rho$ is finally given by Eq. (41).

The most computationally-intensive elements of this approach correspond to steps i and iii, in which functionals of the microscopic interaction potential have to be integrated over the two-particle excluded volume manifold for the chosen mesogen model. Due to the high dimensionality of this manifold, we chose to follow the lead of Ref. Belli et al., 2014 and compute these integrals using stochastic sampling methods, yielding the further advantage of trivial parallelization; the details of this numerical procedure may be found in Ref. Dussi et al., 2015. We further emphasise that this implementation allows for the full evaluation of the angle-dependant ODF $\psi^{(\rho)}_{\rm eq}(\theta)$, rather than its partial projection on a basis of Wigner D-matrices.Tombolato and Ferrarini (2005); Frezza et al. (2014) The current approach is thus both more efficient and not subject to truncation errors.

In the following, we refer to the number of MC integration steps used in both stage i and stage iii by $N_{\rm MC}$, and use a finite grid of fixed size $N_{\theta}=250$ for polar angles $\theta$. The statistical errors of all relevant variables are hence entirely controlled by $N_{\rm MC}$, which can be tuned as required to achieve the desired numerical accuracy.

We remark that this numerical framework is particularly well-suited for the study of lyotropic liquid crystals, as the second-virial coefficient $E^{\infty}$ defined in Eq. (13) is generally independent of particle density. Therefore, the full concentration dependence of the ODF in the Parsons-Lee description may be worked out from Eq. (29) with minimal numerical effort using a single computation of the virial integral (Eq. (13)). Nevertheless, all of our previous discussions remain conceptually valid in the case of thermotropic liquid crystals, and our method may therefore also be applied to the investigation of thermotropic self-assembly by evaluating the temperature dependence of $E^{\infty}$. An important caveat however lies in the limited aspect ratios of most experimental thermotropic mesogens, which may require a more accurate description of inter-particle correlations beyond the decoupling approximation.Xiao and Sheng (2013) We thus largely restrict our focus to lyotropic LC phases in the rest of this paper.

Importantly, it should be noted that solving Eq. (29) in the PL framework at step ii requires a preliminary evaluation of the molecular volume $v_{\rm ref}$ of the reference hard-sphere system in order to work out $\eta_{\rm ref}^{(\rho)}$ from Eq. (22). A number of previous studies obtained good agreement with simulation data for hard convex particles by using the original formulation of the PL theory, which simply amounts to setting $v_{\rm ref}=v$ in Eqs. (22), (28) and (29). However, this assumption has been found to lead to significant underestimations of the osmotic pressure for systems of non-convex particles.Williamson and Jackson (1998) A simple explanation for this discrepancy may be provided by invoking the concept of effective molecular volume,Abascal and Lago (1985) which we introduce below.

Let $\Omega\equiv V\times SO(3)$ be the total configuration space accessible to the particles, with $SO(3)$ the group of rotations in Euclidean 3D space, and $(P,P^{\prime})\in\Omega^{2}$ an arbitrary pair of identical hard particles. The mathematical definition of the excluded volume manifold $\Omega_{\rm exc}$ which may not be penetrated by particle $P^{\prime}$ in any configuration due to the presence of $P$ reads as

[TABLE]

The effective molecular volume $v_{\rm exc}$ of particle $P$ is then defined as

[TABLE]

with $\mu$ the standard Lebesgue measure for subsets of $\mathbb{R}^{3}$. It is easy to show that $v_{\rm exc}$ satisfies $v_{\rm exc}\geq v$, which reflects the smaller volume available to non-convex particles with respect to their convex counterparts at given density $\rho$, and therefore accounts for the higher osmotic pressures obtained in numerical simulations for the former systems. Varga and SzalaiVarga and Szalai (2000) thus proposed to assign this effective molecular volume to the volume of the reference hard-sphere particles,

[TABLE]

which yields the so-called modified Parsons-Lee (MPL) approximation for hard non-convex particles.

A general numerical scheme to compute this effective volume for arbitrary particle models may then be obtained by discretising $V$ into a mesh of size $N_{x}\times N_{y}\times N_{z}$ and evaluating the excluded volume manifold $\Omega_{\rm exc}$ by MC sampling, similarly to the integration procedure outlined above. Assuming the first particle $P$ to be fixed at the origin of the reference frame, this process amounts to generating a number $N$ of uncorrelated configurations for particle $P^{\prime}$ by drawing random center-of-mass positions $\mathbf{r}^{\prime}\in V$ and particle orientations $\mathcal{R}^{\prime}\in SO(3)$, and checking $P$ and $P^{\prime}$ for hard overlaps. In the absence of any such overlaps, all the mesh points contained within the generated configuration for particle $P^{\prime}$ are guaranteed to lie outside of the excluded volume manifold, and may therefore be pruned.

This algorithm thus enables one to efficiently obtain a discretised representation of $\Omega_{\rm exc}$ for a wide range of particle models, with a numerical precision determined by both the chosen mesh resolution and the number of MC steps used. The corresponding reference volume $v_{\rm ref}$ may finally be determined from Eqs. (44) and (45) as the sum of the volumes of all mesh points comprising $\Omega_{\rm exc}$ at the end of the calculation.

All effective excluded volumes used in this paper were computed employing a cuboidal integration box of size $V=L_{x}\times L_{y}\times L_{z}$; $L_{x}$, $L_{y}$ and $L_{z}$ were chosen to yield the minimal volume containing all possible overlapping configurations of two particles, prescribing the center of mass of the reference particle $P$ to the origin and its long axis to $\mathbf{e}_{z}$. The mesh dimensions were then set to $N_{x}\times N_{y}\times N_{z}=150\times 150\times 1500$ and a number $N=1\times 10^{8}$ of random sampling steps was used, yielding a statistical dispersion of less than 1% for all computed values of $v_{\rm ref}$ and ensuring a negligible dependence on varying mesh resolutions.

III Results & discussion

We now illustrate the versatility of our method by determining the nematic and cholesteric phase properties of a variety of mesogenic particles, both chiral and achiral, in the (modified) Parsons-Lee approximation. An interesting first test system to assess its reliability is provided by the coarse-grained helical particle model first introduced in Ref. Frezza et al., 2013, and subsequently studied extensively using both direct computer simulationsFrezza et al. (2013); Kolli et al. (2014a) and different theoretical approaches.Frezza et al. (2014, 2013); Belli et al. (2014); Dussi et al. (2015) Such particles are defined as chains of $N_{s}$ fused hard spheres of radii $\sigma$ placed along a helical backbone of microscopic pitch $p$, radius $r$ and contour length $l$, thus allowing for tunable levels of shape chirality and non-convexity, as illustrated in Fig. 2. Similarly to these previous investigations, we here restrict our study to particles with a right-handed symmetry, and use the shorthand $r_{a}p_{b}l_{c}$ to refer to helices with microscopic radius $r=a$, pitch $p=b$ and contour length $l=c$ in units of $\sigma$.

III.1 The isotropic-cholesteric phase transition

We start by applying our approach to the determination of the isotropic-to-cholesteric phase transition for some of the hard helices above. We first aim to characterise the onset of local nematic order in such systems, and leave the description of the longer-range cholesteric structures arising from particle chirality to Sec. III.3.111We therefore use the words “nematic” and “cholesteric” interchangeably in this section and the next, as our perturbative framework imposes that cholesteric and uniaxial nematic phases be undistinguishable at the local scale. We quantify the degree of local alignment through the usual nematic order parameterde Gennes and Prost (1993)

[TABLE]

The binodal points delimiting the stability ranges of the isotropic and nematic phases can be worked out from the osmotic pressure $\Pi$ and chemical potential $\mu$ of the system, which in our case take the form

[TABLE]

with $\mathscr{F}_{0}$ the undistorted free energy (Eq. (28)), $\Sigma$ the rotational entropy per particle,

[TABLE]

and $b_{2}^{\infty}$ the angle-averaged second virial coefficientLee (1987)

[TABLE]

The binodal coexistence concentrations can then be obtained by equating Eqs. (47) and (48) in the two phases and solving the resulting set of coupled equations via a relaxed multivariate Newton-Raphson scheme,Froberg (1969) using the uniform limit of the ODF for the isotropic phase,

[TABLE]

Some results of our model, specifically the osmotic pressures and nematic order parameters, are summarised in Fig. 2 for different hard helices, including the limiting case of an achiral linear chain of hard spheres (LCHS) for $r=0$. We note that the pressures predicted by our approach are found to be in significantly better agreement with simulation data than those obtained by the theory of Ref. Frezza et al., 2013 in both the cholesteric and isotropic phases for all helices considered. This is because this previous approach assumed the effective molecular volumes $v_{\rm exc}$ of the particles to be given by those of the corresponding LCHS systems with equal contour lengths, for which analytical expressions have been derived.Varga and Szalai (2000) This approximation therefore neglects the contributions of the non-linear molecular backbones to the overall shape non-convexity of the helices, and thus leads to sizeable underestimates of the associated osmotic pressures. These effects may however be accounted for by using the general procedure outlined in Sec. II.4 to work out the full effective molecular volume of each of the chosen particles; the results of both approaches are found to be virtually identical in the case of the LCHS system, for which our numerical excluded volume converges to the analytical expression of Ref. Varga and Szalai, 2000.

Interestingly, we also remark that the quantitative agreement of the MPL theory with simulation results in terms of the nematic order parameter $S$ is somewhat less convincing, and that the inclusion of the full molecular excluded volume does conversely not lead to systematic improvements for the prediction of the I-N coexistence range. This observation seemingly mirrors the findings of Ref. Varga and Szalai, 2000 for systems of long LCHS’s; however, while the MPL treatment was in that case reported to lead to a significant overestimate of the nematic transition density, our results for $r_{.2}p_{2}$ helices rather underestimate this quantity with respect to the simulation data of Ref. Frezza et al., 2013. This would appear to suggest that the influence of particle non-convexity on nematic self-assembly is highly non-trivial, and may be too subtle to be fully captured by the simple hard-sphere analogy of the PL approximation.

III.2 The Frank elastic constants

We now move on to the applications of our approach to the microscopic calculation of the Frank elastic constants, as outlined in Secs. II.3 and II.4. The determination of the LC elastic constants from direct molecular simulations is a notoriously challenging task, owing to the large scale of the associated nematic fluctuations in the weak deformation framework. Therefore, their reliable measurement from the spatial fluctuations of the director field requires the use of very large simulation boxes, which leads to similar technical limitations as mentioned for the direct simulation of cholesteric phases. We thus expect our method to provide a general and efficient way to work out the Frank elastic constants of a wide range of particle models, and to assess their sensitivity to microscopic mesogen structure and thermodynamic state.

However, a consequence of these computational difficulties is the scarcity of the available numerical data to which our results may be compared, as simulation measurements of the elastic constants have so far been limited to a few state points in the phase diagram of highly simplified particle models.Wilson (2005) We therefore only report in Fig. 3 our results in the case of hard spherocylinders with aspect ratios $5\leq x\leq 20$, for which simulationRuzicka and Wensink (2016); Fischermeier et al. (2014); Allen and Frenkel (1988) and theoretical results from both FMTWittmann, Marechal, and Mecke (2015) and DFTSomoza and Tarazona (1989a) have been reported.

We note that our data obeys $K_{3}>K_{1}>K_{2}$ for all particles studied, as expected from general mean-field calculations for calamitic particles.Kröger and Ilg (2007) Furthermore, our results are found to be in reasonable quantitative agreement with the predictions of FMT, although increasing discrepancies are observed in the locations of the I-N transition with decreasing aspect ratios; both approaches are however in comparable agreement with the simulation values of Frenkel and BolhuisBolhuis and Frenkel (1997) for all reported aspect ratios. While both methods provide a good match for the limited simulation data available for the Frank elastic constants, our method predicts a significantly stronger dependence of $K_{1}$ and $K_{2}$ on particle density than estimated by FMT; this deserves to be investigated further, and more extensive simulation studies involving a larger number of state points would therefore be very useful in enabling the assessment of the validity of the different theoretical descriptions.Ruzicka and Wensink (2016)

We also remark that our results are virtually identical to those of Somoza and Tarazona,Somoza and Tarazona (1989a) who made use of a more involved DFT approach based on a weighted-density approximation.Somoza and Tarazona (1989b) This is due to the fact that the full DCF obtained in Ref. Somoza and Tarazona, 1989a reduces to our effective expression Eq. (35) in the functional subspace of Eq. (4), and therefore leads to an equivalent description of nematic director fluctuations in the framework of Sec. II.3. We finally note for the sake of completeness that other results using different variants of DFT have been reported for these systems,Lee (1989); Poniewierski and Hołyst (1990) but are not reproduced here due to their wide dispersion. We further expect the results of FMT, which corresponds to the most sophisticated analytical effort to date in the microscopic description of hard-particle LCs,Mederos, Velasco, and Martínez-Ratón (2014) to provide the most reliable theoretical benchmark currently available for the validation of our approach in this case.

III.3 The cholesteric pitch of hard helices

We finally turn our focus to larger length scales by comparing the cholesteric pitches computed using our method against those reported by previous microscopic theoriesFrezza et al. (2014); Dussi et al. (2015) in order to assess its reliability in predicting the emergence of cholesteric structures at the macroscopic scale. In keeping with the discussions of Sec. III.1, we begin by assessing the dependence of the pitch on the inclusion of particle non-convexity in the MPL approximation. We thus compare in Fig. 4 the results obtained using for $v_{\rm ref}$ both an analytical expressionFrezza et al. (2013) of the real particle volume $v$ and the effective molecular volume $v_{\rm exc}$ computed using the procedure of Sec. II.4 for different hard helices.

We note that the MPL theory does not qualitatively alter the cholesteric behaviour of the systems, and mostly shifts the density dependence of the pitch to lower particle concentrations owing to the smaller volume available to non-convex mesogens. We however emphasise that this effect may lead to significant differences in the handedness and magnitude of the predicted pitches at given particle density, contrary to the suggestion of Ref. Frezza et al., 2014; this is due to the fact that the Parsons-Lee prefactor $G_{\rm PL}$ implicitly affects the equilibrium ODF (Eq. (29)) involved in the calculations of $K_{2}$ and $k_{t}$, even though its contribution to the effective DCF (Eq. (35)) explicitly cancels out at all densities in the ratio of these two quantities (Eq. (41)).

We further remark that the results of Sec. III.1 point to potential deficiencies in the description of the local order parameters — and thus the equilibrium ODFs — worked out by both formulations of the PL approximation for hard helices, which may in turn hinder the reliability of the macroscopic pitches reported by this and previous density functional investigations of these systems; the magnitude of this effect is difficult to assess systematically in the absence of more extensive simulation studies of their nematic and cholesteric behaviour. In order to consistently compare our results with Refs. Belli et al., 2014; Dussi et al., 2015; Frezza et al., 2014, and thus critically assess the effects of our perturbative approach and numerical integration techniques, we choose in the following to restrict ourselves to the canonical form of the Parsons-Lee theory.

Firstly, we compare our results against those predicted by the full-functional approach of Ref. Belli et al., 2014 for weakly chiral systems — corresponding to large equilibrium cholesteric pitches. Fig. 5 summarises our results in terms of the equilibrium cholesteric wavenumber $q_{\rm eq}$ in the case of hard helices with equal microscopic pitch and radius but different contour lengths, for which large equilibrium pitches $\left|\mathcal{P}_{\rm eq}\right|\gtrsim 500\,\sigma$ have been reported.Dussi et al. (2015) We obtain excellent agreement between the two methods throughout the whole cholesteric stability range of all particles studied, providing quantitative evidence of the validity of the perturbative treatment of chirality for these systems. We further note that the computational cost of our method is roughly two orders of magnitude lower than that of our implementation of Ref. Belli et al., 2014, thanks to its use of the full angle-dependant virial coefficient (Eq. (13)) instead of its more expensive, pitch-dependent projections on a finite basis of Legendre polynomials.

We compare in Table 1 the predictions of our approach to those of Refs. Dussi et al., 2015 and Frezza et al., 2014 for a variety of helices at the onset of cholesteric stability. The macroscopic pitches computed by our method are found to be within 15% of those reported by Ref. Dussi et al., 2015 for all particles studied, with the larger relative discrepancies generally corresponding to the shorter pitches, i.e. the more strongly twisted phases. This observation is fully consistent with the limitations of our weak deformation framework, and is further discussed in the next paragraphs.

The correspondence between our results and those of Ref. Frezza et al., 2014 is somewhat more tenuous, especially for larger values of $\mathcal{P}_{\rm eq}$ — a fact all the more puzzling given that both our approaches are based on the same perturbative description of the cholesteric phase. These differences may most likely be attributed to finite-order truncation effects in the Wigner matrix expansions introduced in Refs. Tombolato and Ferrarini, 2005; Frezza et al., 2014 for the computation of $K_{2}$ and $k_{t}$, as small discrepancies in these two coefficients may lead to significant variations of the corresponding equilibrium pitch Eq. (41) in the weak chirality limit. We therefore expect our method, which allows for the calculation of all relevant quantities with a tuneable degree of precision and obviates further numerical approximations, to be more reliable in these cases.

We finally investigate the validity of our perturbative treatment in the limit of strong chirality in Fig. 6 by focusing on systems for which the predicted cholesteric pitches are of the order of a few dozen particle diameters, corresponding to the shortest values reported in Ref. Dussi et al., 2015. Quantitative agreement with the full-functional description of Ref. Belli et al., 2014 is still found to be very satisfactory, which indicates that the effects of the low-order $q$-expansion of the free energy (Eq. (40)) and of the change in local order induced by the cholesteric twist are both of limited consequence for the systems studied.

We nevertheless note that the original framework of Ref. Belli et al., 2014 also relies on the assumption of local phase uniaxiality (Eq. (4)), and therefore still qualifies as a weak-deformation account of cholesteric order. A theoretical route to incorporate the deviations from local cylindrical symmetry induced by macroscopic twist in the approach of Ref. Belli et al., 2014 has subsequently been outlined by Allen,Allen (2016) but has to our knowledge not yet been implemented in numerical studies. These effects are however only expected to be substantial in the limit of short pitches,Priest and Lubensky (1974) and may therefore be of limited practical relevance for the study of most experimentally-realistic systems.

Conversely, it is well known that chiral mesogens are intrinsically biaxial, and that particle interactions at the microscale may thus give rise to an additional level of biaxial organisation in cholesterics beyond this simple phase-induced symmetry breaking.Harris, Kamien, and Lubensky (1997, 1999) In this bottom-up picture, the cholesteric pitch of the system is also expected to be strongly contingent on its degree of local biaxial order, leading to a generally intricate causal relationship between phase chirality and biaxiality.Dhakal and Selinger (2011) Consequently, the use of a uniaxial ansatz for the local description of a cholesteric LC may result in a significant underestimation of its equilibrium pitch. This applies to all studies discussed in this section, and may therefore hamper the reliability of the methods described therein even in the limit of weak deformations. Extending the current density functional framework to explicitly include long-range biaxial correlations would allow for a direct and systematic investigation of the link between particle and phase biaxiality, and of its macroscopic effects on cholesteric behaviour.

III.4 Soft threaded rods: towards more realistic particle models

We finally remark that our method does not rely on any assumptions as to the nature of the microscopic interaction potential $U$ in Eq. (11), and is therefore applicable to any short-range particle potential in the pairwise additive approximation. While the previous sections focused on the case of pure steric repulsion, more experimentally-realistic particle models are likely to involve additional contributions from electrostatic and/or dispersion interactions.Kleman and Lavrentovich (2007) Following Ref. Wensink, 2014, we now consider a system of hard spherocylinders of length $L$ and diameter $D$ decorated with a helical charge distribution of microscopic pitch $p$ located at their surface. To mimic a uniform linear distribution, we discretise the charges in a number $N_{p}$ of identical interaction sites uniformly spread along the helical thread. The full pair interaction energy of two particles $P_{1}$ and $P_{2}$ then reads as

[TABLE]

where the set intersection sign denotes the potential overlap region of the hard spherocylinders, and $u$ represents the electrostatic interaction potential per unit particle length.Wensink (2014) $r_{ij}$ then corresponds to the distance between interaction sites $i\in P_{1}$ and $j\in P_{2}$, which implicitly depends on the center-of-mass separation vector $\mathbf{r}_{12}$ and respective orientations $\mathcal{R}_{1},\mathcal{R}_{2}$ of the two particles. Following Ref. Wensink, 2014, we choose to account for electrostatic repulsion in the Debye-Hückel formalism, introducing a finite cutoff distance $r_{c}$ for computational efficiency,

[TABLE]

with $\lambda_{D}$ the Debye length of the system. In the following, we set the cutoff length to $r_{c}=10\,\lambda_{D}$ and number of patches to $N_{p}=50$, which we verified to yield negligible truncation and discretisation effects for our results, and use $u_{0}=1/D$ for consistency with the values of Ref. Wensink, 2014.

Our numerical approach then enables us to directly investigate the validity of the framework of Ref. Wensink, 2014, in which a similar perturbative description was used, but where a number of approximations were introduced to work out the full pair interaction potential $U$ and Frank parameters $K_{2}$ and $k_{t}$ semi-analytically for this particle model.Wensink and Jackson (2011) For consistency with the approach used therein, we now revert to the original form of the Onsager theory and thus set $G_{\rm PL}=1$ for the rest of this section. The comparison of the equilibrium pitches predicted by our method to those reported in Ref. Wensink, 2014 is then summarised in Fig. 7 for helical charge distributions with a variety of microscopic pitches.

We observe no evidence of qualitative agreement in terms of either magnitude or handedness of the predicted pitches for any of the particles studied, despite the common use of the Onsager-Straley theoretical basis in both our approaches. Furthermore, we find that the inclusion of electrostatic repulsion moves the isotropic-cholesteric transition to lower volume fractions than in the case of purely steric interactions, with the onset of cholesteric order for all systems considered occurring at $\eta\simeq 2\%$ compared to the simulation value of $\eta\simeq 7\%$ obtained for hard spherocylinders with $L/D=50$.Bolhuis and Frenkel (1997) This effect is consistent with the findings of Ref. Stroobants, Lekkerkerker, and Odijk, 1986 for uniformly-charged rods in the Onsager framework, but contrasts with the higher value $\eta\simeq 10\%$ reported in Ref. Wensink, 2014. The origin of these discrepancies most likely lies in some of the simplifying assumptions introduced in Refs. Wensink, 2014; Wensink and Jackson, 2011 to obtain tractable expressions for all microscopic integrals of the form of Eqs. (13) and (32). The magnitude of the resulting disparities illustrates the difficulty of deriving reliable analytical approximations in microscopic descriptions of the cholesteric phase, due to the strong dependence of such systems on the fine structural and thermodynamic properties of their constituent particles.

We also note that a number of simulation results on closely related systems have recently been reported.Ruzicka and Wensink (2016); Honorato-Rios et al. (2016); Kuhnhold and Schilling (2016) However, the first of these studiesRuzicka and Wensink (2016) focused on fairly short spherocylinders in a regime of large Debye lengths, leading to effective particle aspect ratios too low to be accurately described by Onsager-like theories,Stroobants, Lekkerkerker, and Odijk (1986) while the other twoHonorato-Rios et al. (2016); Kuhnhold and Schilling (2016) introduced an explicit dependence of the pair potential $U$ on particle density, which as discussed in Sec. II.4 requires significantly higher computational costs to be probed using our method. We nevertheless remark that Refs. Honorato-Rios et al., 2016; Kuhnhold and Schilling, 2016 also found that the onset of cholesteric organisation was shifted to significantly lower volume fractions by the effects of electrostatic repulsion.

IV Conclusion

We have described a novel density functional method to systematically investigate the link between microscopic structure and macroscopic liquid-crystalline properties in chiral and achiral uniaxial nematic phases. The use of numerical integration to work-out generalised virial coefficients greatly improves its applicability and reliability over that of previous analytical approaches, while the implementation of Straley’s perturbative framework reduces its computational cost to a fraction of that of Ref. Belli et al., 2014 for cholesteric systems. It furthermore circumvents the need for the partial decomposition of the density field in a basis of rotational invariants, which greatly increases its tractability over that of previous numerical descriptions.Tombolato and Ferrarini (2005); Frezza et al. (2014); Belli et al. (2014); Dussi et al. (2015)

The Frank elastic constants computed by our model for hard spherocylinders are found to be consistent with previous numerical and theoretical investigations, while its predictions for the cholesteric pitch of hard helices are in excellent quantitative agreement with those obtained by Ref. Dussi et al., 2015 for all systems studied. The latter confirms that the additional “weak chirality” assumptions introduced by the perturbative approach are not particularly limiting, and still work well down to surprisingly short pitch lengths. We further reiterate that experimental lyotropic cholesterics typically have pitches that are very large compared to the particle diameters, and are therefore expected to be fully amenable to such a treatment.

However, we wish to highlight that there are a number of key approximations in most implementations of the density functional theory of cholesterics that we feel could potentially limit its quantitative accuracy, and need further consideration. Firstly, there are potential deficiencies in the Parsons-Lee rescaling in the case of strongly non-convex mesogens for the prediction of the onset of nematic order, which are not obviously improved upon by modifying the theory to account for the real effective molecular volume of the particles. This likely points to the need for a more sophisticated description of inter-particle correlations in these cases, following the considerations of Sec. II.2. We note that an approach involving the explicit third-order virial contribution to the excess free energy was also implemented with some success for such systems, albeit under stringent symmetry assumptions.Kolli et al. (2014b); Cinacchi et al. (2016) It would be interesting to investigate whether its generalisation to configurations of freely-rotating particles may lead to significant improvements over the MPL approximation in our case.

Secondly, there is the assumption of local uniaxiality, and we wish to further emphasise the potential importance of long-range biaxial correlations in systems of chiral mesogens. The explicit angle dependence of our excluded volume integrals provides a natural framework for the introduction of phase biaxiality, which may shed some light on the systematic underestimations of the cholesteric twist obtained by density functional methods in previous studies.Varga and Jackson (2006); Dussi and Dijkstra (2016) The explicit account of biaxial order will also represent a first step towards the description of more complex chiral nematic phases beyond the cholesteric symmetry, such as the so-called twist-bendDozov (2001) and screw-likeKolli et al. (2014a) nematic phases. However, recent studies suggest that the accurate treatment of the former may require a microscopic theory of nematic elasticities beyond the perturbative Poniewierski-Stecki framework,De Gregorio et al. (2016) while the coupling between orientational and translational degrees of freedom of the latter would call for further modifications of our symmetry assumptions.Kolli et al. (2014a) Thus, in our opinion, significant efforts are still required to achieve a unified field-theoretical description of the chiral nematic state.

Hence, we stress that numerical simulations of simple cholesteric model systems along the lines of Refs. Dussi and Dijkstra, 2016; Ruzicka and Wensink, 2016; Kuhnhold and Schilling, 2016 are of particular importance to allow the quantitative validity of these different approximations to be probed. However, the determination of cholesteric pitches in these previous studies relies on rather involved numerical procedures, and their degree of mutual consistency remains largely unascertained due to their focus on different particle models and parameter ranges. We therefore hope that this paper will stimulate further investigations of the respective reliabilities of these different methods, and generate a reference set of data suitable for the thorough assessment of density functional descriptions of the cholesteric phase.

We finally point out that the use of Monte-Carlo techniques for the computation of all microscopic integrals further allows for the simple treatment of non-rigid particle models, e.g. by drawing mesogen configurations from an ensemble of representative structures at every integration step as part of the random sampling process. This will enable us to directly investigate the influence of particle flexibility on phase chirality, an issue so far largely unexplored despite its essential experimental relevance.Cherstvy (2008); Dogic and Fraden (2006)

Acknowledgements.

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 641839. The authors would like to acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work. http://dx.doi.org/10.5281/zenodo.22558.