A variational mean-field theory for spatially modulated liquid crystal phases

We generalize a recent Maier–Saupe theory for thermotropic bent-core mesogens²⁴ to determine the phase behavior of curved spherocylinders with diameter D, length L, and radius of curvature R corresponding to a central angle Ψ = L/R (Fig. 1d). We describe a curved spherocylinder with centre-of-mass position R = (X, Y, Z) and orientation Ω = (α, β, γ) as a rigid chain of M segments of length L/M. We find that M = 10 segments is sufficient to account for the particle shape for the range of Ψ that we considered. Each segment i, with centre-of-mass position r_i and orientation $$ {\widehat{\mathbf{u}}}_i$$, is assumed to align preferentially along the local nematic director $$ \widehat{\mathbf{n}}\left({\mathbf{r}}_i\right)$$ via an effective mean-field potential βU(R, Ω). In addition, we employ McMillan’s extension²⁵ to account for possible density modulations along the (global) nematic director. The resulting effective mean-field potential reads

As a variational ansatz for $$ \widehat{\mathbf{n}}\left(\mathbf{r}\right)$$ we employ the nematic director field $$ {\widehat{\mathbf{n}}}_{\mathrm{TSB}}\left(z\mid {\theta }_{\mathrm{a}},{\theta }_{\mathrm{b}},q\right)$$ of a generic N_TSB phase with conical angle θ_a and θ_b, wavenumber q, and the helical axis along the z − direction

We note that the N_TSB phase reduces to an N_TB phase with a circular precession cone when θ_a = θ_b = θ₀, whereas for θ_a = θ_b = π/2 the circular cone reduces to a flat circle resulting into an N^* phase as the precession of the nematic director reduces to a simple twist around the phase director. If either θ_a or θ_b vanishes, the elliptical cone collapses onto a flat isosceles triangle and an N_SB phase is obtained with nematic director field $$ {\widehat{\mathbf{n}}}_{\mathrm{TSB}}\left(z\mid {\theta }_0,0,q\right)=\sin \left({\theta }_0\right)\sin \left(qz\right)\mathbf{j}+\sqrt{1-{\sin }^2\left({\theta }_0\right){\sin }^2\left(qz\right)}\mathbf{k}$$ which differs from the expression $$ {\widehat{\mathbf{n}}}_{\mathrm{SB}}\left(z\right)=\sin \left({\theta }_0\sin \left(qz\right)\right)\mathbf{j}+\cos \left({\theta }_0\sin \left(qz\right)\right)\mathbf{k}$$ given by Dozov¹⁴. However, for small conical angles the two expressions for the nematic director fields are indistinguishable, see Supplementary Note 1. Finally, for θ_a = θ_b = 0 the cone becomes a simple line, and the N_TSB simply reduces to a uniaxial N phase. Hence, all the above-mentioned nematic phases are limiting cases of the N_TSB phase, see Fig. 2. In the following, we distinguish N_TB, N_TSB, and N_SB states based on their ellipticity e = θ_b/θ_a: we label N_TB the states with e > 0.8, N_TSB the states with 0.2 < e < 0.8, and N_SB the states with e < 0.2. We note that the thresholds e = 0.2 and e = 0.8 — necessary because of the statistical error on e — are arbitrary, and every state with ellipticity e ≠ 0, 1 is, in principle, an N_TSB state.

As βU(R, Ω) only depends on the Z-component of R with period p, we restrict all integrations over R to integrations over Z ∈ [0, p]. The onset of orientational and/or positional order corresponds to a change of entropy per particle $$ \mathrm{\Delta }\mathcal{S}/N=-k_{\mathrm{B}}\int \mathrm{d}Z\mathrm{d}\mathrm{\Omega }f\left(Z,\mathrm{\Omega }\right)\log \left[16{\pi }^{3}f\left(Z,\mathrm{\Omega }\right)/q\right]=⟨U\left(Z,\mathrm{\Omega }\right)⟩/T-k_{\mathrm{B}}\log \left(16{\pi }^3/qQ\right)$$and to a change of energy per particle ΔU/N = 〈U(Z, Ω)〉/2. Using Eq. (2), we obtain the change of free energy per particle relative to the isotropic state

Minimizing ΔF with respect to S, τ, λ, and the variational parameters θ_a, θ_b, and q of the nematic director field $$ {\widehat{\mathbf{n}}}_{\mathrm{TSB}}\left(z\mid {\theta }_{\mathrm{a}},{\theta }_{\mathrm{b}},q\right)$$ yields the equilibrium phase at temperature k_BT/ϵ. We note that since the nematic director $$ {\widehat{\mathbf{n}}}_{\mathrm{TSB}}\left(z\mid {\theta }_{\mathrm{a}},{\theta }_{\mathrm{b}},q\right)$$ already imposes a periodicity with pitch length p = 2π/q in the system, p must be a multiple of the periodicity λ of the density modulations, i.e.p = nλ with $$ n\in \mathbb{N}$$. As the numerical minimization of ΔF always yields non-integer values of n close to 2 in spatially modulated phases, we impose n = 2 in all considered cases.

In Fig. 3a we show the resulting phase diagram from Maier–Saupe theory for a system of curved spherocylinders as a function of Ψ and inverse temperature βϵ, where we set α = 0.05. At low curvature, i.e. small Ψ, the I phase transforms into a uniaxial N phase and subsequently into an N_TB phase upon increasing βϵ. However, the stability range of the N phase shrinks with increasing particle curvature, eventually disappearing at Ψ ≳ 1.2. Our results show that deformations of the nematic director field become more pronounced with increasing particle curvature Ψ until the I-N phase transition is replaced by a direct I-N_TB transition^20,24,26.

Fig. 3

Phase diagram of curved spherocylinders.

Phase behaviour of hard curved spherocylinders as (a) predicted by theory, as a function of inverse temperature βϵ and particle curvature Ψ, and (b) obtained from simulations, as a function of packing fraction η and Ψ. Both phase diagrams exhibit an isotropic (I) phase at low βϵ and η, respectively. For small particle curvatures Ψ the I phase transitions into an uniaxial nematic (N) phase upon increasing βϵ or η. The stability range of the N phase decreases with Ψ and disappears at Ψ ≳ 1.2. Upon increasing βϵ or η further, the twist-bend nematic (N_TB) phase transforms into a twist-splay-bend smectic (Sm_TSB) phase, and eventually into a splay-bend smectic (Sm_SB) phase. The splay deformations in the nematic director field are accompanied by density modulations. Source data are provided as a Source Data file.

Moreover, our generalized Maier–Saupe theory predicts that upon increasing βϵ further, the nematic director field exhibits not only twist and bend modulations but also splay deformations, resulting in a twist-splay-bend phase. Upon lowering the temperature further, the twist deformations are gradually replaced by splay deformations, eventually yielding a splay-bend phase at sufficiently high βϵ. In particular, we find two distinct regions of splay-bend phases, one at low particle curvature and one at high curvature, the latter transforming into a re-entrant twist-splay-bend phase with increasing βϵ. Intriguingly, our theory predicts that the onset of splay deformations is accompanied by density modulations. We therefore refer to these phases as twist-splay-bend smectic (Sm_TSB) and splay-bend smectic (Sm_SB) phases rather than N_TSB and N_SB phases. This finding is also supported by the observation that in the case of α = 0, i.e. without McMillan’s extension to describe density modulations, the phase diagram displays only stable I, N, and N_TB phases. The N_SB phase is thus unstable in the case of a homogeneous density field. Furthermore, we find that varying the value of α > 0 results only in a temperature shift of the transition from N_TB to the Sm_TSB and Sm_SB phase, whereas the amplitude of the density modulations remains unaffected. This confirms that splay deformations of the nematic director field are inherently associated with the appearance of density modulations, as already speculated by De Gennes and Meyer, four decades ago^27,28.

Monte Carlo simulations of hard curved spherocylinders in bulk and sedimentation

To test the predictions of our Maier–Saupe theory, we study the bulk phase behavior of hard curved spherocylinders with L/D = 19 and varying particle curvature Ψ using NPT-MC simulations, i.e. the number of particles N, pressure P and temperature T are fixed. The resulting phase diagram is shown in Fig. 3b as a function of Ψ and packing fraction η. The phase behavior of this athermal lyotropic system is driven by η that plays a similar role as the inverse temperature βϵ in our Maier–Saupe theory for thermotropic systems. Using this analogy, the comparison of the topology of the theoretical with the computational phase diagram is remarkable. Our simulations confirm the I-N-N_TB phase sequence at low particle curvature that transforms into Sm_TSB and Sm_SB phases upon increasing η, as well as a direct I-N_TB transition at high particle curvature transforming into Sm_TSB and Sm_SB phases with increasing density. Our simulations reveal that the splay modulations are accompanied by density modulations in agreement with our predictions from Maier–Saupe theory. In the Supplementary Note 4, we present a mapping of the theoretical and computational phase diagrams using the dependence of the global nematic order parameter S_g on temperature and packing fraction, and we provide a comparison of the orientational order parameters, pitch and conical angles as a function of thermodynamic state, showing quantitative agreement between the predictions of the Maier–Saupe theory and simulations. We note that the main qualitative difference is the absence of the re-entrant Sm_TSB phase in the simulated phase diagram. It is important to mention here that due to hysteresis effects and slow equilibration it is impossible to accurately determine the regions of stability and the first- or second-order nature of the phase transitions of the N_TB, Sm_TSB, and Sm_SB phases in simulations of hard curved spherocylinders. However, the simulations consistently show that the phase transformation from an N_TB to an Sm_SB phase proceeds via an intermediate Sm_TSB phase as twist deformations are gradually replaced by splay deformations of the nematic director field. For example, in Fig. 4 we show typical configurations of an Sm_SB phase, an Sm_TSB phase, and an N_TB phase along an expansion of a system of hard curved spherocylinders of length L/D = 19 and opening angle Ψ = 1.31 from packing fraction η = 0.406 to packing fraction η = 0.367. To fully characterise the phases in Fig. 4, we plot the scalar order parameter S(z), the nematic director field $$ \widehat{\mathbf{n}}\left(z\right)$$, and the density field ρ(z) in Supplementary Fig. 12, 13 and 14. In addition, we also plot the polarisation vector $$ \widehat{\mathbf{m}}\left(z\right)$$ along with the bend distortions in the director field described by the bend vector $$ \widehat{\mathbf{b}}\left(z\right)=\widehat{\mathbf{n}}\times \left(\nabla \times \widehat{\mathbf{n}}\right)/\parallel \widehat{\mathbf{n}}\times \left(\nabla \times \widehat{\mathbf{n}}\right)\parallel $$. We clearly observe that $$ \widehat{\mathbf{m}}\left(z\right)$$ is always anti-parallel to $$ \widehat{\mathbf{b}}\left(z\right)$$, demonstrating that these modulated phases are driven by bend deformations coupled to polar ordering spontaneously arising from the packing constraints of banana-shaped particles.

Fig. 4

Splay-bend smectic Sm_SB, twist-splay-bend smectic Sm_TSB, and twist-bend nematic N_TB phases.

Typical configurations of (a) an Sm_SB phase with conical angles θ_a ~ 0.76 and θ_b ~ 0 at packing fraction η = 0.394 (box size 55.5D × 30.3D × 47.7D), b an Sm_TSB phase with conical angles θ_a ~ 0.87 and θ_b ~ 0.51 at packing fraction η = 0.371 (box size 52.0D × 33.4D × 49.1D), and (c) an N_TB phase with conical angles θ_a ~ θ_b ~ 0.83 at packing fraction η = 0.354 (box size 49.9D. 3D × 49.3D) along an expansion of a system of hard curved spherocylinders of length L/D = 19 and opening angle Ψ = 1.31.

Finally, we perform simulations on a system of hard curved spherocylinders with L/D = 19 and Ψ = 0.99 under gravity with a gravitational length l_g = k_BT/mg = 7.5D parallel to $$ \widehat{\mathbf{z}}$$ with a hard wall at the bottom at z = 0, m denoting the buoyancy mass of the rods and g the gravitational acceleration. The resulting configuration, presented in Fig. 5a, shows the full I-N-N_TB-Sm_TSB-Sm_SB phase sequence in a single system. In particular, we observe a continuous transition from an N_TB phase with a director field precessing on a circular cone at the top of the sediment, via an Sm_TSB phase where the precession cone becomes elliptic, towards an Sm_SB phase where the elliptical cone reduces to a flat triangle at the bottom of the sample. Hence, the transition from the N_TB to the Sm_SB phase occurs via a continuous flattening of the precession cone of the nematic director field from a right circular cone to an isosceles triangle via a continuous range of elliptical precession cones corresponding to a range of Sm_TSB phases. In the Supplementary Note 6, we compare the equation of state obtained by integrating the density profile of the sediment with the equation of state obtained from bulk simulations, showing good agreement.

Fig. 5

A system of curved rods in gravity.

a Typical configuration from a simulation on a system of hard curved spherocylinders with an aspect ratio L/D = 19 and central angle Ψ = 0.99 in a gravitational field along $$ -\widehat{\mathbf{z}}$$ with a gravitational length l_g = 7.5D and a hard wall at z = 0. The phase sequence I-N-N_TB-Sm_TSB-Sm_SB is observed from the top to the bottom in the sediment. b The x-, y-, and z-components of the nematic director field $$ \widehat{\mathbf{n}}\left(y\right)=\left(n_x\left(y\right),n_y\left(y\right),n_z\left(y\right)\right)$$ in subsequent slabs of the system of thickness 5D, showing a continuous transition from splay-bend deformations to twist-bend deformations with increasing altitude z, corresponding to lower pressure and density. The measured values of the nematic director field are denoted by the varying symbols, whereas the lines denote a fit of the simulation data using Eq. (3). The precession cones correspond to the fits of the nematic director field, showing a continuous transition from an N_TB phase with a director field precessing around a circular cone at the top of the sediment, via an Sm_TSB phase where the precession cone becomes elliptic, towards an Sm_SB phase where the elliptical cone reduces to a flat triangle at the bottom of the sample. Source data for (b) are provided as a Source Data file.

Rationalising the topology of modulated liquid crystal phases

Our extensive simulations show that our generalized Maier–Saupe theory effectively predicts, despite its simplicity, the stability of twist-bend, twist-splay-bend, and splay-bend phases of curved spherocylinders. Our microscopic theory is solely based on the tendency of particles to align their particle shape to the local nematic director field $$ \widehat{\mathbf{n}}\left(\mathbf{r}\right)$$. To rationalize our findings we calculate the integral curves of the nematic director field $$ \widehat{\mathbf{n}}\left(z\right)$$, i.e. curves r(z) to which the nematic director $$ \widehat{\mathbf{n}}\left(z\right)$$ at z is tangent, such that $$ {\mathbf{r}}^{\prime }\left(z\right)=\nu \widehat{\mathbf{n}}\left(z\right)$$ with ν a proportionality constant. We thus find

For a generic twist-splay-bend (TSB) phase the integral in Eq. (5) cannot be evaluated analytically. However, the tendency of particles to align their profiles to the nematic director field at low temperatures or high densities, corresponds to a tendency to match their particle curvature with the curvature of its integral curves. Given a curve r(z) with tangent $$ {\mathbf{r}}^{\prime }\left(z\right)=\nu \widehat{\mathbf{n}}\left(z\right)$$ of constant norm $$ \parallel {\mathbf{r}}^{\prime }\left(z\right)\parallel =\nu $$, we can reparametrize it by its arc length as r(s = νz). In this reparametrization the tangent to the curve $$ {\mathbf{r}}^{\prime }\left(s\right)=\left(\partial \mathbf{r}\left(z\right)/\partial z\right)z^{\prime }\left(s\right)=\widehat{\mathbf{n}}\left(s\right)$$ has unit norm, and we can calculate the curvature of the integral curve as $$ \kappa \left(s\right)=\parallel {\widehat{\mathbf{n}}}^{\prime }\left(s\right)\parallel $$. Going back to the original parametrization, we obtain $$ \kappa \left(z\right)=\parallel {\widehat{\mathbf{n}}}^{\prime }\left(z\right)\parallel /\nu $$. Using this formula, we obtain the following curvature for the integral curve of the nematic director field of a generic TSB phase,

In the presence of splay deformations, i.e. θ_a ≠ θ_b, κ_TSB(z) is a non-trivial periodic function of z. In Fig. 6a–b, we show the curvature κ_SB(z) of the integral curves of the nematic director field of various splay-bend (SB) phases with θ_a = θ₀ and θ_b = 0 from theory and simulations, along with their density profiles ρ(z), measuring the probability of finding a particle at z. The periodicity of κ_SB(z) agrees well with that of $$ \log \rho \left(z\right)$$, i.e. minus the effective mean-field potential felt by the particles. If the curvature of the integral curves of the nematic director field matches that of the particles, the particles can optimally align their shape to the local nematic director field, resulting into a lower potential energy or higher free volume and thus a higher local density ρ(z), whereas geometric frustration arises when the curvature of the integral curve deviates from that of the particles, leading to a higher potential energy or lower free volume and hence a lower local density. We thus conclude that the presence of periodic density modulations in TSB and SB phases, and hence Sm_TSB and Sm_SB phases, can be explained geometrically by the coupling of particle curvature to the non-uniform bend deformations of the director field, but also by the coupling of density to splay deformations^27,28.

Fig. 6

Rationalization of the spatial modulations of the density and the nematic director fields.

a–b Modulations of the curvature κ(z) of the integral curves of the nematic director field (top), and of the logarithm of the density profile ρ(z) (bottom) corresponding to minus the effective potential acting on the particles, in splay-bend smectic states (Sm_SB) predicted by our theory (a) and found in Monte Carlo simulations (b) for a system of hard curved spherocylinders of various curvatures Ψ and inverse temperatures βϵ or packing fractions η as labelled. c Pitch length p versus the conical angle θ₀ of various twist-bend nematic (N_TB) phases of hard curved spherocylinders of various curvatures Ψ ∈ [0.5, 2] from theory (top) and simulations (bottom), along with the relationship of Eq. (7) (pink line). Source data are provided as a Source Data file.

On the other hand, in the case of an N_TB phase the translational symmetry of the elastic deformations is preserved, see Supplementary Note 2, and we expect the curvature of the integral curves of the nematic director field to be uniform. Intriguingly, for N_TB phases with θ_a = θ_b = θ₀ the integration of Eq. (5) can be carried out explicitly, yielding the following expression $$ {\mathbf{r}}_{\mathrm{TB}}\left(z\right)={\mathbf{r}}_{\mathrm{TB}}\left(z_0\right)+\frac{\nu }{q}\sin {\theta }_0\left[\sin \left(qz\right)-\sin \left(q{z}_0\right)\right]\mathbf{i}-\frac{\nu }{q}\sin {\theta }_0\left[\cos \left(qz\right)-\cos \left(q{z}_0\right)\right]\mathbf{j}+\nu \cos {\theta }_0\left[z-z_0\right]\mathbf{k}$$ for the integral curve of the nematic director field, i.e. a helix of period $$ p\nu \cos {\theta }_0$$. Hence, we can impose that the integral curve has the same pitch p as the N_TB phase — or, equivalently, that $$ \mathbf{r}\left(z\right)\cdot \widehat{z}=z$$ — by setting $$ \nu =1/\cos {\theta }_0$$. As expected, the curvature of Eq. (6) reduces to a uniform value $$ {\kappa }_{\mathrm{TB}}\left(z\right)=\left(q/\nu \right)\sin \left({\theta }_0\right)=q\sin \left(2{\theta }_0\right)/2$$ independent of z. Using the conjecture that particles tend to match their curvature with the one of the integral curves of the nematic director field, we impose κ_TB = 1/R, and obtain the simple relation

Methods

Variational Maier–Saupe theory

To solve our generalized Maier–Saupe theory, we minimize the free energy in Eq. (4) with respect to the variational parameters of the director- and density-field. Calculating the free-energy difference per particle βΔF/N is trivial except for the partition function Q, i.e. an integral of the form

with $$ f\left(Z,\mathrm{\Omega }\right)=\exp \left[-\beta U\left(Z,\mathrm{\Omega }\right)\right]$$. We assume that the function f(Z, Ω) is a periodic function in Z with period p. To evaluate the integrals of the form (8), we first transform the original space of integration to a 4-dimensional hypercube [−1, 1]⁴ using the transformation $$ \alpha =\pi \left(\xi +1\right),\cos \beta =\eta ,\gamma =\pi \left(\mu +1\right),Z=\frac{p}{2}\left(\psi +1\right)$$ with Jacobian ∣J∣ = pπ²/2. The integral in Eq. (8) becomes

which we solve using an $$ \mathcal{N}$$-points Gauss-Legendre quadrature. If {ξ_i}, {η_i}, {μ_i}, and {ψ_i} are the $$ \mathcal{N}$$ node points in [−1, 1], i.e. the $$ \mathcal{N}$$ roots of the $$ \mathcal{N}$$-th Legendre polynomial, and {w_ξ,i}, {w_η,i}, {w_μ,i}, and {w_ψ,i} are the associated weights, we can approximate the integral I of Eq. (9) by

Using the approximation of Eq. (10) and $$ f\left(Z\left(\psi \right),\mathrm{\Omega }\left(\xi ,\eta ,\mu \right)\right)=\exp \left[-\beta U\left(Z\left(\psi \right),\mathrm{\Omega }\left(\xi ,\eta ,\mu \right)\right)\right]$$, we evaluate the partition function in Eq. (4) using an $$ \mathcal{N}$$ = 32 points Gauss-Legendre integration.

Once Q is calculated, we minimize βΔF/N at given βϵ in the 6-dimensional space of parameters S, τ, n, θ_a, θ_b and q = 2π/p by means of a Covariance Matrix Adaptation Evolution Strategy (CMA-ES) using the Python library in https://pypi.org/project/cma/.

Monte Carlo simulations

We study the phase behavior of a system consisting of hard curved spherocylinders using Monte Carlo simulations in the NPT ensemble, i.e. at fixed number of particles N, pressure P and temperature T. We employ an orthorhombic simulation box of sides L_x, L_y, and L_z and apply periodic boundary conditions. We perform a sequence of MC cycles consisting of N + 1 MC moves. Each MC move consists of a particle move with probability ~N/(N + 1), and a volume move with probability ~1/(N + 1). In a particle move, a random roto-translation of a randomly picked particle is proposed and accepted if it does not generate overlaps with other particles. In a volume move, a random variation of a random side of the simulation box is proposed, and the system is compressed or expanded accordingly. If the compression/expansion does not generate overlaps between the particles, the move is accepted with a probability

where $$ \mathrm{\Delta }V=V^{\prime }-V$$ denotes the change in volume. We perform simulations on a system of N = 2048 hard curved spherocylinders, and initialize all simulations from a nematic configuration. We measure a wide set of observables during the simulation, like density, uniaxial nematic order parameters, smectic order parameters, etc. When the system reaches equilibrium as monitored from the observables, we characterize the system’s configuration. We find that 10⁸ MC cycles are typically sufficient for equilibration.

To study a system of hard curved spherocylinders in a gravitational field, we perform Monte Carlo simulations in an NVT ensemble, i.e. we fix the number of particles N = 8192, volume V, and temperature T. We implement a hard wall at the bottom at z = 0 and apply periodic boundary conditions in the x − and y − direction. Each MC cycle consists of N attempts to randomly rotate and translate a randomly selected particle. If a particle roto-translation leads to an overlap with one of the particles or with the hard wall, the move is rejected, otherwise it is accepted with a probability

with $$ \mathrm{\Delta }{U}_{\mathrm{g}}=\left(z^{\prime }-z\right)/l_{\mathrm{g}}$$ the change in potential energy due to gravity, z and $$ z^{\prime }$$ the old and new z-coordinate of the displaced particle, and l_g the gravitational length of the system.