In silico study of liquid crystalline phases formed by bent-shaped molecules with excluded-volume type interactions

We have numerically studied a liquid composed of achiral, bent-shaped molecules built of tangent spheres. The system is known to spontaneously break mirror symmetry, as it forms a macroscopically chiral, twist-bend nematic phase [Phys. Rev. Lett. 115, 147801 (2015)]. Here, we have examined the full phase diagram of such liquid and observed several phases characterized by orientational and/or translational ordering of molecules. Apart from conventional nematic, smectic A, and the above-mentioned twist-bend nematic phase, we have identified antiferroelectric smectic A phase. For large densities and a high degree of molecule's structural bend, another smectic phase emerged, where the polarization vector rotates within a single smectic layer. These results were confirmed using both Monte Carlo and molecular dynamics simulations.


Introduction
Comprehension of spontaneous mirror symmetry breaking (SMSB) phenomenon is essential for understanding fundamental aspects such as life itself [1][2][3] . Investigations revolving around the origins and corollary of SMSB stand simultaneously at the forefront and frontiers of contemporary soft matter science 4,5 . Recently, SMSB was observed in liquid crystals (LCs), where supramolecular chiral and polar structures are assembled from achiral, bent-shaped (curved) molecules 6 . Bent-shaped molecules (BSMs) tend to form, i.a., smectic and nematic phases, whereas the latter have recently attracted a resurgence of worldwide interest due to novel findings 7,8 . One of them is the discovery of twist-bend nematic (N TB ) phase [9][10][11] (Fig. 1), which constitutes the first example of SMSB in a liquid state with no support from long-range spatial ordering. The characteristic feature of N TB is the emergence of heliconical (1D modulated) structure of nanoscale pitch with a ground-state exhibiting a degenerate sign of chirality (ambidextrous chirality), which means that each helix handedness is equiprobable. Further advances in the studies on N TB subsequently sprouted into a broad sub-field of liquid crystal research dedicated to novel nematic phases of non-standard symmetries [12][13][14] .
It is noteworthy to point out that remarkable features of nematic phases formed by BSMs are not limited solely to the phenomenon of SMSB. The list is vast but just to name a few 17,18 : negative bend-splay elastic anisotropy (i.e. K 33 − K 11 < 0), giant flexoelectricity, high Kerr constants, and large rotational viscosity. The BSMs are also said to be the perfect candidates for optically a Institute of Theoretical Physics, Jagiellonian University, Łojasiewicza 11, 30-  biaxial fluid, namely biaxial nematic (N B ), predicted by Freiser in 1970 19 . It is so because the introduction of the bend deformation in the molecular structure naturally specifies a secondary direction for orientational ordering, which is perpendicular to the long molecular axis 20,21 .
All this makes LCs formed by BSMs, per se, a riveting field to explore both from experimental and theoretical points of view 15,16,18,[22][23][24][25][26] . In particular, substantial efforts were undertaken to decipher how the molecule's structure influences the emergence of individual mesophases 21,27,28 . One of valuable tools in such studies are computer simulations, such as Monte Carlo sampling (MC) and Molecular Dynamics (MD), which have a long-standing tradition in the investigation of structural properties of LC phases using a diverse range of models ranging from simplified lattice systems to fully atomistic ones [29][30][31] . Generally, the anisotropic shape of BSMs can be approximated in the simulations either by V-like or C-like bent objects (reflecting the C 2v symmetry) composed of i.a., needles 32-34 , spheres [35][36][37][38] , or spherocylinders [39][40][41][42] . It has been shown that those systems, each one separately, produce a wealth of interesting nematic and smectic phases. However, especially appealing are the systems where the supramolecular chirality emerges out of pure entropy-driven interactions.
Herein, we have extended the work of Greco and Ferrarini 38 where it was shown that for a model system of BSMs (C-shapelike), composed of eleven tangent spheres -see Fig. 2, it is possible to observe the stabilization of N TB arising from the packing entropy in MD simulations. We aim to provide an in-depth analysis of this system by altering the BSMs' curvature and packing densities in order to obtain the full phase diagram for the liquid state * . To achieve it we performed both MC and MD simulations.
The article is organized as follows. In Sec. 2 the system under study is described together with the details concerning the MC and MD simulations used to investigate its properties. The results of the simulations are given in Sec. 3 where they are compared with the current literature. In Sec. 4 we provide conclusions drawn from the acquired data.

The model system
BSMs were modeled by eleven tangent spheres of diameter σ placed on a circular arc with a variable bend angle (see Fig. 2). For χ = 180 • , the molecule reduces to a linear chain of adjoint spheres , while for χ = 0 • it forms a semi-circle. To study lyotropic phase transitions in a system composed of such molecules, we examined two types of excluded-volume interaction. The first one is hard-core interaction 44 , while the second is Weeks-Chandler-Anderson (WCA) 45,46 soft repulsion of each pair of spheres, de- Fig. 2 Left: Molecular structure of N TB forming mesogen CB7CB 9 rendered via QuteMol 43 . Right: Model bent-shaped molecule with C 2v symmetry composed of eleven tangent spheres whose centers are equidistantly distributed on the arc. The bend angle χ is defined as the angle formed by two lines, tangent to the arc at the centers of end-point spheres. Unit vectorsâ andb define the molecular frame, whereâ equates with molecular long axis andb is parallel to the molecule's two-fold (C 2 ) symmetry axis.
fined as where r mn is the distance between the interaction centers and ε is the energy scale. The free parameters of the model σ , ε can be conveniently incorporated into the reduced pressure P * = Pσ 3 /ε and reduced temperature T * = k B T /ε. For hard-core repulsion, the equation of state depends only a single parameter -the P * /T * = Pσ 3 /(k B T ) ratio 44,47 . For WCA interaction the system is sensitive to P * changes, while the dependence on T * is weak 48 .
The emergent phases were studied as a function of bending angle χ and packing density η. The latter one is a natural parameter for lyotropic systems and enables us to compare the two interaction models 48 . For both of them, the main contribution to Helmholtz free energy is the packing entropy: F ≈ −T S (the formula is exact for hard-core repulsion), so the phases with the highest entropy minimize F and thus are promoted.

Simulation methods
The hard-core interactions were simulated using Monte Carlo sampling in NPT ensemble 49 . The system was built of N = 3000 molecules. A single MC cycle consisted of N molecule moves and a volume move. To perform a molecule move, a single molecule was sampled at random. Then, the molecule was translated by a random vector and rotated by a random angle around a random axis. The move was accepted provided it did not introduce an intersection. For volume moves all three lengths of an orthorhombic simulation box with periodic boundary conditions and the positions of molecules were scaled by three independent logarith-mically distributed factors. Moves introducing overlaps were rejected and the rest of them was accepted according to Metropolis criterion. The extents of random moves were selected such that the acceptance ratio was around 15%. The initial configuration was a diluted antiferroelectric sc crystal, which was then compressed or expanded to a target density. The equilibration was very sluggish due to concavities of molecules and required up to 2×10 8 of full MC cycles. The production runs to sample ensemble averages consisted of 3×10 7 cycles. To accelerate the simulations, we parallelized them using domain division 50 .
The MD simulations of N = 2025 molecules, interacting via the WCA potential in the NPT ensemble at various values of the pressure P * in the interval 0.2 ≤ P * ≤ 2.0, were performed using the open-source software LAMMPS 51 . All relevant simulation aspects were set up according to the framework provided in 52 . Equilibration and production runs were conducted for 10 7 timesteps each.

Order parameters
Distinctive features of a particular phase can be characterized by appropriate order parameters and a corresponding degree of ordering. In nematic and smectic phases, the molecules' long axesâ tend to orient, on average, along a preferred direction called directorn. Such kind of alignment can be quantified by the secondrank order parameter 53,54 where the sum is over all N molecules in the system, P 2 is the Legendre polynomial of degree 2,â i is a molecular axis vector of i-th molecule (cf. Fig. 2) and . . . denotes the time average over uncorrelated system snapshots. If P 2 = 0 orientations are isotropic, P 2 = 1 is a perfect order, and for P 2 = −1/2 all molecules are perpendicular ton. Alternatively, P 2 can be computed by diagonalising the second rank Q-tensor 55 whereâ i,α andâ i,β are the Cartesian components ofâ i and δ αβ is the Kronecker delta. The nematic order parameter P 2 (for a single snapshot) is associated with the largest-modulus of a nondegenerate eigenvalue of Q andn is the corresponding eigenvector. In this formulation, zero-energy diffusive motion ofn does not have effect on P 2 value.
The density modulation is measured by smectic order parameter 53,54 where ρ is number density, k is a wavevector of density undulation and r i are molecules' centers of mass. Typically, for non-tilted smectic phases k n. For τ = 0 there is no density modulation, whereas for τ = 1 the molecules form perfect layers with a spacing equal to 2π/|k|.
Total polarization can be computed as a sum of molecule po-larization vectors As it will be seen from the results, the non-zero polarization can be present when one restricts the summation to a plane perpendicular ton, however, at the same time, it can vary in the direction ofn and give zero net polarization when summed over the whole system. To account for such variations, one can calculate the norm M(C j ) = M(C j ± ∆C/2) for molecules in a narrow shell between two planes r ·n = C j ± ∆C/2 and sum them over the whole system Here, N b = L/∆C is the number of shells and L is the system length in the direction ofn.

Phase overview
The system exhibits a variety of liquid phases: isotropic liquid (Iso), nematic (N), twist-bend nematic (N TB ), non-polar smectic A (SmA), antiferroelectric smectic A (SmAP A ) and an unidentified polar smectic (SmX). The phase diagram for hard-core interactions is shown in Fig. 3. WCA soft repulsion recreates all the phases with nearly identical phase borders. All order parameters described in the previous section are gathered in Fig. 4 as a function of the bend angle χ and the packing fraction η. Finally, system snapshots for all identified phases are depicted in Fig. 5. For each value of χ and sufficiently low η, the molecules form ordinary isotropic liquid as shown by low values of all order parameters indicating lack of any long-range correlations [see Fig. 5(a)]. Around η = 0.45 the system crystallizes into several types of hexagonal close-packed crystals, which is beyond the scope of this work. The boundary of the Iso phase moves upwards in η for a decreasing opening angle χ. It can be explained by the fact, that a lower bend angle results in less prolate molecules, for which the entropic gain from orientational ordering is lower.

Nematic phases
The ordinary nematic phase N can be observed for χ > 110 • , as indicated by a sharp jump of P 2 value from near zero to over 0.5 on the phase boundary [see also Fig. 5(b)]. A growing η results in an increase of P 2 to almost 0.9 for χ ≈ 180 • , which is significantly higher than a typical range [0.3, 0.7] for nematics observed in experiments 56 . In this range of parameters (χ, η), no density modulation is observed. Moreover, two smallest eigenvalues of Q are comparable, meaning that the nematic is uniaxial, which is further supported by a vanishing net polarization M ≈ 0 (result not shown). Layer polarization m is also near zero, so the phase is homogeneous.
The situation changes for a phase forming over the Iso phase for χ < 110 • and over the N phase in the range χ ∈ [110 • , 140 • ], for η values between 0.29 and 0.35 [see Fig. 5(c)]. No density modulation suggests that is is a nematic, however with P 2 ∈ [0.3, 0.5], lower than for the ordinary N phase in the system. The net polarization is zero, while m lies in the range [0.4, 0.8], indicating high polar ordering in planes perpendicular ton. A closer inspection reveals that this is twist-bend nematic phase N TB , observed in many studies of bent-core molecules 38,42,57,58 (see also Fig. 1). In this phase, bothn and polarization vary in space in a heliconical pattern. Assuming the global director is aligned with the z axis, the director field may be parameterized aŝ The spatial dependence of M is similar, with M ⊥ n. N TB phase has two parameters -the conical angle θ , which is the tilt of the director from the z axis, and a pitch p defining the periodicity of the structure. The parameters vary in phase spacep ∈ [25σ , 45σ ] and θ ∈ [25 • , 35 • ]. With increasing packing fraction η or decreasing bend angle χ the period p shortens and θ becomes more ob-tuse. N TB can be predicted using phenomenological lowest-order Oseen-Zocher-Frank free energy 59-61 where subsequent terms represent, respectively, splay, twist and bend deformations of director field and K ii (i = 1, 2, 3) are Frank elastic constants. All three Frank elastic constants have to obey Ericksen inequalities 62 , i.e.: K ii ≥ 0 (i = 1, 2, 3), otherwise perfect alignment (homogeneous nematic) would not correspond to a state of minimum energy. However, if one were to consider scenario where K 33 < 0, then it would occur that bend deformation is energetically favorable and reference state, i.e. homogeneous nematic, becomes unstable 63 . As there does not exist a structure with constant bend without topological defects, is has to be combined with at least one of: splay or twist deformations. If K 11 > 2K 22 , twist is preferred over splay 63 . As pointed by Shamid et al. 58 , the effective K 33 constant can become negative below a certain critical temperature, when the coupling of polarization with bend deformation is introduced, supporting the assumption that N TB is inherently polar. The phase can be also understood via geometrical reasoning 15,16,26,64 . For curved spherocylinders, one can predict a correct relation between molecule bend angle, pitch p and conical angle θ matching the spherocylinder curvature with an integral curve of a director fieldn 42 .

Smectic phases
For η > 0.34, smectic phases are formed. The N phase is adjacent to an ordinary SmA phaseτ is over 0.5 and strata are clearly visible [see Fig. 5(d)]. Further compression narrows the spread of molecule mass centers in the direction normal to the layer and τ approaches unity for a high η. Typically, P 2 value experiences a jump on the N-SmA boundary from around 0.4-0.6 to over 0.8 57 , however in some systems it is small or absent 40 . Here, no sudden change of P 2 is observed. Molecule polarizations are chaotic within the layers, as indicated by a low value of m and they are orthogonal ton.
A distinct smectic phase is formed over N TB phase for χ > 100   reveals that adjacent layers have an opposite (antiferroelectric) polarization, allowing us to classify the phase as SmAP A . The onset of polar order may be attributed to the entropic gain of translational motion within a layer, however, it does not explain the antiferroelectric arrangement of layer polarization. Lansac et al. 40 attributed such a configuration to the more compatible vibrational motion of two molecules from adjacent layers when they have antiparallel directions compared to parallel ones, resulting in once again in a higher entropy of the former system. SmAP A -SmA boundary moves to higher densities for less bent molecules -lower curvatures discourage polar order, so a stronger compression is required to achieve it. In particular, for χ > 160 • , a crystallization occurs before the formation of SmAP A † . Note that N TB always forms a polar smectic when compressed. It can be expected, since N TB is already polar. A further insight into the connection between SmAP A and N TB can be given by the analysis ofn(x, y, z) spatial dependence. The numerical data is well described by a director field n(x, y, z) =    sin(θ ) cos(2πz/p) 0 thus the emergent SmAP A phase can also be described as splaybend smectic (Sm SB ). There, the director is constant in xy plane and oscillates within xz plane aroundẑ direction, coinciding with the direction k of density modulation. Althoughn leans from k within a layer, which would suggest type C smectic, averaging over the whole layer givesn k indicating that this is SmA, as previously stated. The director fields (7), (9) are special cases of a more general director field modulation with all of splay, twist, and bend deformationŝ where the angles θ a and θ b may be different. Fixing θ a = θ and continuously changing θ b from 0 to θ interpolates between director fields for N TB and Sm SB . This intermediate splay-twistbent smectic Sm STB phase has been observed for bent spherocylinders 42 . In our model, however, it is absent. † Notably, the crystal is polar for χ ∈ [160 • , 170 • ] and non-polar above χ = 170 • .
Highly bent molecules form yet another polar smectic phase, as indicated by high τ and m values, however of an unknown type, denoted as SmX [see Fig. 5(f)]. Despite trying both compression and expansion from various initial configurations (dense crystal, diluted crystal, type C smectic, lower-density equilibrium system, and others) we were not able to obtain a structure without irregular domain walls. As it can be seen from Fig. 5, the polarization within each stratus rotates with y coordinate and topological defects are visible. The same effects are present for both hard and soft interactions. Increasing the system size, which allows a full rotation of the molecule polarization vector has not eliminated the defects and in neither case they form a regular pattern as in e.g. blue phases 65 . However, the blue-phase-like structures cannot be completely disregarded since the lattice constant of the defects can be of the order of tens or even hundreds of molecule lengths, and the system size allowing such long periods is far beyond the limit of our computational resources. Another candidate is anticlinic, antiferroelectric type C smectic (SmC A P A ) 66 . There, both the director tilt and polarization alternate between layers. Remarkably, we have observed its formation for a very acute bend angle χ < 70 • and a high density η > 0.4 (see Fig. 6).
BSMs with highly acute bend angle are still scarce and not fully understood 66,67 . Additionally, chemical data shows that increasing the degree of bend diminishes the LC phase stability to the point where the molecular bend is so acute that the compound becomes non-mesogenic 68 . Therefore, χ < 90 • falls outside of the scope of this work. Nevertheless, further investigation of the unidentified SmX phase can be a potential goal for future studies.

Conclusions
We have systematically studied the phase diagram of the liquids formed by bent-shaped molecules composed of tangential spheres with centers placed on an arc. The molecules were characterized according to adjustable bend angle χ, and the system of them was considered under different packing densities η. In (χ, η) parameters space, we have identified isotropic, nematic, twist-bend nematic, smectic A, and antiferroelectric smectic A phases. For the highest density and the largest χ we have observed a smectic phase, denoted as SmX, forming polar layers with in-layer rotating polarization. We were, however, not able to avoid defects and domain walls. For unrealistically large bend, SmX phase transformed to anticlinic, antiferroelectric type C smectic, where both the director tilt and polarization alternate between layers. In comparison to a previous study for bent spherocylinders 42 , we have not observed splay-twist-bend smectic. The obtained results were the same for both: MC with hard-core interactions and MD using short-range WCA repulsion.

Conflicts of interest
There are no conflicts of interest to declare.