Jagiellonian University 


Jagiellonian University Institute of Physics 
Department of Statistical Physics 

Doctoral thesis 

Lattice Models of Biaxial, Tetrahedratic 
and Chiral Order 

Karol Trojanowski 
Supervised by prof. dr hab. Lech Longa 
Cracow, Summer 2013 

Peer-reviewed papers 
• 	
Karol Trojanowski, David W. Allender, Lech Longa and Łukasz Kuśmierz. (2011). Theory of Phase Transitions of a Biaxial Nematogen in an External Field. Molecular Crystals and Liquid Crystals, 540(1), 59-68. 

• 	
Karol Trojanowski, Grzegorz Pająk, Lech Longa and Thomas Wydro. (2012). Tetrahedratic mesophases, chiral order, and helical domains induced by quadrupolar and octupolar interactions. Physical Review E, 86(1), 011704. 

• 	
Lech Longa and Karol Trojanowski. (2013). Ambidextrous Chiral Domains in Nonchiral Liquid-crystalline Materials. Acta Physica Polonica B, 44(5), 1201-1207. 

• 	
Karol Trojanowskia and Lech Longa. (2013). Synchronization of Kuramoto Oscillators with Distance-dependent Delay. Acta Physica Polonica B, 44(5), 991. 



Conference talks and posters 
• 	
March 2009. 35th Conference of the Middle European Cooperation in Statistical Physics, Pont-a-Mousson, France. Poster “Synchronization of Phase-coupled Oscillators With Distance-dependent Delay”. 

• 	
July 2010. 23rd International Liquid Crystal Conference, Kraków, Poland. Poster “Biaxial phase of interacting quadrupoles in an external field”. 

• 	
February 2011. 11th European Conference on Liquid Crystals, Maribor, Slovenia. Oral presentation “Biaxial phase in nematic liquids under influence of external fields”. 

• 	
September 2011. 24th Marian Smoluchowski Symposium on Statistical Physics, Zakopane, Poland. Poster “Effects of external fields in the interacting quadrupoles model of biaxial nematics”. 

• 	
November 2011. 12th International NTZ-Workshop on New Developments in Computational Physics, Leipzig, Germany. Oral presentation “Frequency-locking and pattern formation in the Kuramoto model with Manhattan delay”. 

• 	
September 2012. 25th Marian Smoluchowski Symposium on Statistical Physics, Kraków, Poland. Poster “Lattice simulations of tetrahedratic and chiral order”. 


• 	
March 2013. 38th Conference of the Middle European Cooperation in Statistical Physics, Trieste, Italy. Poster “Tetrahedratic mesophases, chiral order, and helical domains induced by quadrupolar and octupolar interactions”. 

• 	
August 2013. 26thMarianSmoluchowskiSymposiumonStatisticalPhysics,Kraków, Poland. Oral presentation “Frequency-locking and pattern formation in the Ku-ramoto model with Manhattan delay”. 

• 	
September 2013. 13th European Conference on Liquid Crystals, Rodos, Greece. Oral presentation “Chiral structures of nonchiral liquid-crystalline materials”. 



Abstract 

The developments over recent years in the research of liquid crystals with non-standard symmetry fuel intensified endeavors, as scientists turn their attention to exotic liquid crystalcompounds, suchasbent-coremesogens,ferrocenemesogensand,recently, flexible dimers, which, accordingtomanyreports, producephasesunlikeanythingstudiedbefore. Among the most notable properties of these phases is ambidextrous chirality – where in the nematic phase the sample forms domains of opposite handedness. Furthermore, the biaxial nematic phase continues to garner interest. This phase had been long missing, but has been supposedly discovered in a bent-core liquid crystal in 2004 – a finding which was subsequently criticized, but spurred a heated dispute over the experimental feasibility of stabilizing the biaxial nematic phase, which has also lead to discussion whether the phase could be field-stabilized (and if this was actually the case in the famous experiments). Further results hint also at long-range tetrahedratic ordering and the curious phenomenon of spontaneous chiral symmetry breaking and forming of chiral twisted macroscopic states. The twist-bend nematic phase, which has been recently discovered for flexible dimers, appears to exhibit helical-conical twists of pitch at the nano-scale, apparently spanning only several molecules, which is particularly surprising, since in most of the chiral phases, such as the cholesteric phase, the pitch is hundreds to thousands of molecules long. Since a satisfying theoretical account of these effects in the form of a microscopic model has not been yet given, the present thesis aims to recreate these phenomena in the framework of a generalized lattice dispersion model, studied by Monte Carlo simulations and by application of mean-field theory. 
After a general introduction to liquid crystals and their symmetries in Chapter 1, in Chapter 2 the model is defined through the formalism of symmetrized irreducible spher<BR>ical tensors, which we use as symmetry-adapted building blocks. In the dispersion inter-action potential, which can be represented as coupling of molecular multipolar moments of two molecules, we consider the molecular quadrupolar moment Q(ˆ
Ω), which consists of a cylindrically-symmetric uniaxial part and a biaxial part of rectangular-cuboidal symmetry, and the molecular octupolar moment T(3), which is a spherical tensor of the 
2 
symmetry of a regular tetrahedron. Two intermolecular coupling terms are constructed out of the respective multipolar moments involving coupling of moments of the same order exclusively. Additionally, a term quadratic in a generalized external vector field is introduced, where the molecular diamagnetic or dielectric polarizability tensor is assumed proportional to the molecular quadrupolar moment. Finally, two terms which account for coupling of the quadrupolar and octupolar moments to the intermolecular vectors are proposed. The full Hamiltonian considers intermolecular interaction between nearest neighbors on a simple cubic lattice with periodic boundary conditions. This general model is subsequently studied in three limit cases, the model of a biaxial nematic in an external field in Chapter 3, the model with quadrupolar and tetrahedratic octupo<BR>lar coupling terms in Chapter 4 and its extension to intermolecular vector coupling in Chapter 5. 
The first considered case, which is the model of a biaxial nematic in an external field, can be regarded as an extension of a well-known model, which was studied extensively by Luckhurst and Romano, who built up from the pioneering work of Freiser in 1970 and Straley in 1974. This model is known to produce the Maier-Saupe transition from the isotropic to the uniaxial nematic phase, as well as a low-temperature second-order transition to the biaxial nematic phase. In the parametrization which we use, a scalar parameter λ mixes two symmetrized irreducible tensors, which correspond to uniaxial and biaxial symmetries to compose the molecular quadrupolar moment Q(ˆa second
Ω) – rank symmetric and traceless tensor. Q(ˆ
Ω) can be regarded as the anisotropic part of the molecular (dielectric or diamagnetic) polarizability tensor and varying λ represents modulating the shape of the molecule. The spectrum of λ covers shapes from uniaxial prolate (rodlike) to uniaxial oblate (discotic) molecules, with biaxial shapes of either tendencyinbetween. AsymmetrytransformationoftheLuckhurst-Romanomodelexists between the prolate and oblate sides of the range of possible values of λ, with a pivot point at a special value, called the self-dual point (at which the model is symmetric with itself), which also corresponds to completely (maximally) biaxial molecules (ones which cannot be unequivocally classified as either prolate or oblate). By an addition of a term quadratic in field and linear in Q(ˆ
Ω) to this model, we study the field-induced effects on the bulk sample of mesogens which exhibit interactions of uniaxial and biaxial symmetry. 
Firstly, the effects of arbitrary non-zero field are predicted through considering the minimum of the field interaction potential in all phases, for both cases of anisotropy of the molecular polarizability (positive or negative) and with respect to varying λ. This analysis reveals a duality transformation of the model which involves switching the sign of molecular anisotropy and the shape tendency (shifting λ about the self-dual point) simultaneously, thus linking the effects encountered for prolate molecules with positive molecular anisotropy to those encountered for oblate molecules with negative molecular anisotropy and, conversely, those encountered for prolate molecules with negative anisotropy to effects present for oblate molecules with positive molecular anisotropy. In the last two cases the field induces biaxiality in the nematic phase. The field effects at the self-dual point need to be treated separately because of the high symmetry involved. 
Subsequently,thephasediagramsofthemodelarestudiedinthespaceofthetemperature and a field parameter, which is the product of the squared field magnitude and the molecular anisotropy. (This allows us to consider the cases of positive and negative anisotropy and varying field magnitude on the same phase diagram by using only one parameter.) The isotropic to uniaxial nematic phase transitions terminates at a critical point for prolate molecules of positive molecular anisotropy and for oblate molecules of negative molecular anisotropy, as the field is increased. For the remaining cases the transition is made continuous at a tricritical point for a sufficiently large magnitude of the field. Using a mean-field high-temperature expansion the critical and tricritical points are studied. Values of critical field and critical temperature are found for both cases in the entire molecular shape spectrum, parametrized by λ. The main result states that when the self-dual point is approached from either side (the molecular biaxiality is increased), the critical field for both the critical and tricritical points is lowered and is zero at the self-dual point. The critical temperature for the both cases is found to increase with λ. Mean-field calculations and Monte Carlo simulations are carried out for selected values of λ and at the self-dual point to study the exemplary phase diagrams for the model in detail. As expected on the grounds of the duality transformation, the shapes of the phase diagrams for prolate and oblate molecules are mutually inverted with respect to the field parameter. In both cases the temperature of the isotropic to uniaxial nematic phase transition is increased with field. The biaxial phase and the transition temperature is not altered by the increased field in the cases of prolate molecules of positive molecular anisotropy and for oblate molecules of negative molecular anisotropy, while for the remaining cases the biaxial phase is also induced in the uniaxial nematic phase and the zero-field transition to the biaxial nematic phase is removed. For maximally biaxial molecules (for λ at the self-dual point) the external field works to reduce the temperature of the transition to the biaxial phase (for either case of molecular anisotropy), which is the only case of field-lowering of the transition temperature for this model. 
The lattice dispersion model with quadrupolar and tetrahedratic octupolar coupling terms was originally postulated by Longa, Pająk and Wydro [1] and our aim here is to extend the research by providing detailed analysis in the terms of temperature scans of order parameters, fluctuation of energy and order parameter susceptibilities for some of the cases which have not been yet considered. In Chapter 4 we study the model for representative values of λ (including the self-dual point) and for two cases of the coupling constant which scales the tetrahedratic interaction term, τ . The first case corresponds to a phase diagram dominated by nematic phases, already obtained in [1], but for which no temperature dependence of the order parameters and susceptibilities has been previously given. The second case leads to a phase diagram with a multicritical point at which six phases: isotropic, nematic uniaxial oblate, tetrahedratic, nematic tetrahedratic prolate, nematic tetrahedratic oblate and chiral nematic tetrahedratic meet. The results for the nematic-dominated phase diagram corroborate the previous result from [1]. The biaxial nematic phase exists only in the vicinity of the self-dual point, while away from it there exist sequences of isotropic, nematic uniaxial, nematic tetrahedratic and chiral nematic tetrahedraticphases,withdecreasingtemperature. Thephasetransitionsassociatedwith spontaneous breaking of tetrahedratic symmetry are manifestly first-order and the phase transitions to the biaxial nematic and chiral phases are second order. For the second case, for which the multicritical point is found, there exists a direct first-order phase transition from the isotropic phase to the chiral nematic tetrahedratic phase at the selfdual point. For the remaining cases the expected purely tetrahedratic phase is not found, as it exists in a too narrow temperature range. Instead, a sequence of phases is found of isotropic, nematic tetrahedratic and chiral nematic tetrahedratic phases, with lowering temperature. As before, the isotropic to nematic tetrahedratic phase transition is first order, while the nematic tetrahedratic to chiral nematic tetrahedratic phase transition is second-order. The chiral nematic tetrahedratic phase is not only chiral, but also spatially homogeneous and biaxial. This might come across as counter-intuitive, as chirality is known to lead to inhomogeneous ground states. However, this model does not explore all of the couplings which can be constructed with the tensors present in the theory. 
In particular, an antisymmetric tensor can be constructed out of the tetrahedratic octupolar molecular tensor and the uniaxial and biaxial parts of the molecular quadrupolar tensor. This leads to the interaction terms involving coupling of the quadrupolar and octupolar tensors to intermolecular lattice vectors, considered in Chapter 5. The quadrupolar version of this interaction is inserted into the model with quadrupolar and tetrahedratic octupolar coupling. The resulting model is given a preliminary investigation, consisting mainly of simulation results in the form of configuration snapshots and commentary. Firstly, a basic discussion of the ground state of the interaction potential leads to a conclusion that the minimum is achieved when two neighboring particles are rotated with respect to each other. This is demonstrated by simulation on a onedimensional chain. We follow with a discussion on how the ground state in two and more dimensions leads to frustration, which results in the formation of complex structures. Simulations on a two-dimensional lattice show that there exists a phase transition from the disordered to the chiral phase (this is the only phase transition in two dimensions, as other degrees of freedom are continuous), in which there appears to be some kind of reordering and two distinct types of structures are present in the lower and higher temperature regimes. One is identified as a bent-wavefront-cholesteric structure, in which cholesteric-like layers bend at right angles, but form no defects. The second is a simpler layeredstructure,akintoatwo-dimensionalcholesteric. Inthreedimensionswehavealso identified two distinct structures in the chiral phase, out of which the low-temperature one is an unidentified complex supramolecular structure, while in the high-temperature structure the molecules form helical twists in a layer-like structure with a pitch spanning only several molecules. In this case the bulk cross-sections bear an uncanny resemblance to the freeze-fracture patterns only recently obtained in experiments on the twist-bend nematic phase. 
In the final chapter we comment on the entirety of the obtained results and on the many cases which have not been investigated in this thesis. 

Acknowledgements 

The Monte Carlo simulation software used in this work is partly inspired by the software originally written by Tomasz Wydro. 
This work is supported by the International PhD Projects Programme of the Foundation for Polish Science within the European Regional Development Fund of the European Union, agreement no. MPD/2009/6. 

The research was carried out with the supercomputer “Deszno” purchased thanks to the financial support of the European Regional Development Fund in the framework of the Polish Innovation Economy Operational Program (contract no. POIG.02.01.00-12023/08). 
x 



Contents 
Papers and conference appearances iii 
Abstract v 
Acknowledgements x 
1 Introduction 1 
1.1 Liquid crystals, typical phases and symmetries . . . . . . . . . . . . . . . 1 
1.1.1 Typical thermotropic mesophases and symmetries . . . . . . . . . . 2 
1.1.2 Molecularchirality........................... 3 

1.2 Non-standardphases.............................. 4 

1.3 Introduction to order parameters in liquid crystals . . . . . . . . . . . . . 6 
1.4 Dispersionforces ................................ 8 

1.5 Historicalresults ................................ 9 

1.5.1 Maier-Saupe model of thermotropic nematics (lattice case) . . . . . 9 
1.5.2 Lebwohl-Lasherlatticemodel..................... 11 

1.6 Alignmenttensor................................ 11 

1.7 Mean-fieldtheoryforlatticemodels ...................... 13 

1.8 MonteCarlosimulations............................ 15 

1.8.1 Metropolisalgorithm.......................... 15 

1.8.2 Rotationaldegreesoffreedom..................... 17 

1.8.3 Parallel sampling and pseudorandom number generation . . . . . . 19 
1.9 Purposeandplanofthesis........................... 20 

2 Generalized dispersion model for bent-core (and related) systems 23 
2.1 Irreducibletensorsandsymmetrization.................... 23 

2.1.1 Symmetrized irreducible tensors for selected symmetry groups . . . 25 
2.2 Definitionofthemodel............................. 26 

3 Dispersion model of biaxial nematics in an external field 31 
3.1 Introduction................................... 31 

3.2 Dispersionmodelsofbiaxialnematics . . . . . . . . . . . . . . . . . . . . 31 
3.2.1 Dispersion model of interacting quadrupoles . . . . . . . . . . . . . 32 
3.2.1.1 Properties of Q(ˆ34
Ω) ...................... 

3.2.1.2 Dualitytransformation ................... 36 

3.2.1.3 Phasesandphasetransitions . . . . . . . . . . . . . . . . 36 
3.3 Dispersionmodelwithexternalfield ..................... 38 

xi 

3.3.1 Model .................................. 39 

3.3.2 Field-inducedeffects .......................... 39 

3.3.2.1 Uniaxiality-promoting effects . . . . . . . . . . . . . . . . 40 
3.3.2.2 Biaxiality-promoting effects . . . . . . . . . . . . . . . . . 41 
3.3.2.3 Field effects for λ <BR>= <BR><BR>
................... 42 

6 
1√ 
3.3.3 Dualitytransformation......................... 42 

3.4 Toolsandmethods............................... 43 

3.4.1 Orderparameters............................ 43 

3.4.1.1 Alignmenttensor....................... 43 

3.4.1.2 Definitionoforderparameters . . . . . . . . . . . . . . . 45 
3.4.1.3 Invariant parameter w <BR>.................... 45 

3.4.2 Mean-fieldmethod........................... 46 

3.4.2.1 Self-consistentapproach................... 47 

3.4.2.2 High-temperature expansion . . . . . . . . . . . . . . . . 48 
3.4.3 MonteCarlosimulations........................ 48 

3.4.3.1 Specific heat and susceptibilities . . . . . . . . . . . . . . 50 
3.5 ResultsfortheLuckhurst-Romanomodel . . . . . . . . . . . . . . . . . . 52 
3.6 Studyingthemodelwithexternalfield . . . . . . . . . . . . . . . . . . . . 58 
3.6.1 Criticalandtricriticalpoints ..................... 58 

3.6.1.1 Mean-field calculation of the critical point . . . . . . . . . 59 
3.6.1.2 Mean-field calculation of the tricritical point . . . . . . . 61 
3.6.1.3 Conclusions and comments on simulations and experiment 63 
3.6.2 Results for λ <BR>=0.3 <BR>........................... 64 

3.6.3 Results for λ <BR>=0.5 <BR>........................... 72 

3.6.4 Results for λ <BR>= <BR><BR>
1√
........................... 74 

6 
3.7 Summary .................................... 78 

3.7.1 Phasediagramsandorderparameters . . . . . . . . . . . . . . . . 78 
3.7.1.1 λ <BR>=0.3 <BR>............................ 78 

3.7.1.2 λ <BR>=0.5 <BR>............................ 79 

3.7.1.3 <BR><BR>
λ = 

1√
............................ 

6 
79 

3.7.1.4 Mean-field study of critical and tricritical points . . . . . 80 
4 Tetrahedraticorderandchiralsymmetrybreakinginadispersionmodel of liquid crystals 81 
4.1 Introduction................................... 81 

4.2 Dispersion model with tetrahedratic coupling (infinite pitch limit) . . . . . 82 
4.2.1 Orderparameters............................ 84 

4.2.2 Phasesofmodel(4.1) ......................... 85 

4.2.3 Scopeofpresentedresearch ...................... 87 

4.3 MonteCarlosimulations............................ 88 

4.3.1 Susceptibilities ............................. 89 

4.3.2 Correlationfunctions.......................... 89 

4.3.3 Determining transition temperatures . . . . . . . . . . . . . . . . . 90 
4.4 Results for τ <BR>=1 <BR>................................ 91 

28
4.5 Results for τ = 15 ............................... 95 

4.6 Discussionandsummary............................ 99 

5 Spatially modulated structures 101 
5.1 Intermolecularvectorcoupling .........................101 

5.2 Consequences of intermolecular vector coupling . . . . . . . . . . . . . . . 102 
5.2.1 One-dimensionalchain.........................104 

5.3 Two-dimensionallattice ............................105 

5.3.1 Simulationsandresults ........................105 

5.4 Three-dimensionallattice ...........................108 

5.5 Discussion....................................110 

6 Summary and conclusions 111 
Bibliography 115 
Chapter 1 

Introduction 
1.1 Liquid crystals, typical phases and symmetries 
Liquid crystals are intermediate states of matter between liquids and crystalline solids, possessing some of the properties of both. Anisotropic molecules, called mesogens, through mutual interaction produce macroscopic states which exhibit ordinary rheology like that of liquids, and macroscopic anisotropic properties observed in crystals, such as anisotropy of the dielectric polarizability, birefringence and Bragg peaks in X-ray and neutron scattering [2]. The liquid properties are due to the retained translational free-dom, while anisotropy is an effect of long-range orientational ordering of the mesogens, which on average align along macroscopic axes of anisotropy. Liquid crystalline states, calledmesophases,occurinsystemsdividedintotwogroups. Inthermotropicliquidcrystals phase transitions between mesophases occur with changing temperature (or altering the pressure). Several phases of different symmetry can be observed in sequence with decreasing temperature for the same compound, before it crystallizes to a solid. On the other hand, there are lyotropic liquid crystals, where the transitions occur with changing concentration of a component. They are, in general, mixtures of different compounds distributed in a solvent and the macroscopic ordering can be produced by alignment of individual molecules, but also of larger structures, such as anisotropic micellae [3, 4]. In the presented work we deal with thermotropic mesophases exclusively. 
The initial discovery of a mesophase is attributed to botanist Friedrich Reinitzer, who in 1888 noticed that the synthesized by himself cholesteryl benzoate exhibits two thermotropic phase transitions in the liquid phase [5]. Subsequently, he performed several experiments involving polarized light, which revealed the anisotropic nature of the newly found state of matter. Puzzled by these discoveries, he contacted Otto Lehmann, who 
1 

initially asserted that the intermediate cloudy phase is in fact crystalline, but later concluded that it is in fact an anisotropic liquid and continued Reinitzer’s research [6, 7]. Because the new phase bore the properties of both liquids and crystalline solids, the term “liquid crystal” was coined and is used to this day, generally in reference to compounds which produce mesophases. Liquid crystal research was not prominent for nearly half a century, but experienced an outburst after World War II, when in the oncoming decades followed many discoveries of new phases and successful theoretical descriptions [2–4, 8– 10] along with numerous practical and commercial applications, such as ubiquitous liquid crystal displays [11, 12] and electro-optical devices [13]. 
1.1.1 Typical thermotropic mesophases and symmetries 
Thermotropic liquid crystal systems above the highest melting point (called the clearing point) are ordinary liquids, characterized by the isotropic symmetry of the O(3) group of proper and improper rotations in 3D space. Lowering the temperature leads to a phase transition to a mesophase, in which directions of anisotropy are established by spontaneous breakdown of rotational symmetry. The simplest example of mesophase, although not the first to be discovered, is the nematic phase. Calamitic (rod-like or cigar-shaped) mesogens with large longitudinal anisotropy (e.g. MBBA or PAA [2, c. 1.2]) possess (on average) an effective axis of cylindrical symmetry, while the intrinsic degrees of free-dom and anisotropy of secondary axes average out due to thermal fluctuations and the molecules align with their principal molecular axes in parallel. Among classical nematogens there are no known cases of compounds which produce stable polar nematic phases [14], hence the average alignment axis, called the director, has no preferred orientation. The (uniaxial) nematic phase has global symmetry of the D∞h point subgroup of O(3), which consists of continuous rotation about a fixed axis and one reflection plane perpendicular to the axis of symmetry. In this case, the global symmetry axis is the director. A phase of identical symmetry is also formed by discotic nematogens (e.g. hexasubstituted triphenylene [15]), which in average resemble flat disks. In this case the molecular axis of symmetry is perpendicular to the face of the disk, so that the macroscopic alignment of faces produces a uniaxial nematic phase. The two varieties are distinguished as uniaxial nematic prolate NU+, for the alignment of long molecular axes, and uniaxial nematic oblate NU−, for the alignment of faces. 
Traditionally mentioned are also layered mesophases called smectics, in which there also occurs partial breaking of the translational symmetry and the distribution of centers of mass exhibits one-dimensional periodicity. Out of the vast variety of smectic phases, we name a few. The simplest smectic, known as smectic A (SmA), shows nematic order 

�Barrett Research Group, Creative Commons Attribution-Share Alike 3.0 Unported under fair use.) 
within periodic layers with the director normal to the layer. In smectic C (SmC) the director is tilted with respect to the layer normal by an angle. 

1.1.2 Molecular chirality 
The phenomenon of molecular chirality in mesogens leads to many new phases. Chiral molecules are not symmetric with respect to a rotoreflection, just like one’s right hand cannot be rotated to match the left and vice-versa. Such molecules are optically active, 
i.e. there is an asymmetry in the interaction of the molecule with the two helicities of light and the interaction can alter the light’s polarization [16]. Homogeneous macro<BR>scopic chirality occurs when a compound, for various reasons, does not possess an isomer of opposite chirality or when the specific synthetization process favors one of them. The other case is that molecules of opposite handedness are equally abundant in a mixture. In this case, domains of opposite chirality can be formed. In terms of microscopic alignment, molecular chirality leads to rotated ground states of two molecules with the same chirality, with a preferred helicity for two left-handed and two right-handed molecules. In homochiral samples or domains, this induces mesoscopic twisted states. Particularly, in chiral nematogens, in a homochiral domain there can exist a macroscopic wave vector along which the perpendicular director rotates, forming a helix with a helicity reflecting the handedness of the molecules. Cross sections perpendicular to the wave vector exhibit nematic ordering, while being rotated with respect to each other. This case describes the cholesteric phase, which was first found in cholesteryls (making it the first liquid crystal phase discovered). 
In some compounds the ground state is frustrated such that the energy minima with respect to different directions cannot be achieved simultaneously and the frustration is relaxed by forming superstructures known as blue phases [17, 18]. Many chiral mesogens also form chiral layered phases, the simplest of which is the chiral smectic (SmC∗) phase, 

�Barrett Research Group, Creative Commons Attribution-Share Alike 3.0 Unported under fair use.) 
inwhichthetilteddirectorissubsequentlyrotatedbetweenneighboringlayers(effectively being constrained to the surface of a cone.) 


1.2 Non-standard phases 
While the uniaxial nematic has one axis of symmetry breaking (one Goldstone mode), in a biaxial nematic the D∞h symmetry is further spontaneously broken and a secondary axis is formed. Because two distinct axes naturally distinguish a third one, normal to the plane of the first two, the resulting phase is expected to have anisotropic properties in three distinct directions. The symmetry of the biaxial nematic phase NB is of the point group D2h, the symmetry of a rectangular cuboid (Fig. 1.4a), which consists of reflections along three perpendicular mirror planes and three twofold symmetry axes perpendicular to the mirror planes. 
By analogy, since uniaxial molecules produce the uniaxial nematic phase, the biaxial phase should be formed by biaxial (brick-like or ellipsoloidal) molecules which align their primary and secondary axes in a long-range fashion. This was notably observed for anisotropic micellae in a lyotropic system by Yu and Saupe [19] and more recently in polymeric systems by Severing et al. [20], but the discovery of a thermotropic mesogen whichproducesthebiaxialnematicphasehasprovedfarmoredifficult. Intriguingly,most nematogens are intrinsically biaxial [21], but produce mostly stable uniaxial phases. The existence of the biaxial nematic phase has been subject to intensive studies [22–31]. 
Fairly new accounts report observation of the biaxial nematic phase in systems of bentcore mesogens [27, 29], as predicted by Teixeira et al. [32], and tetrapode molecules [28]. The discoveries have been notably questioned [33] and alternative explanations of the results have been proposed [34], thus in the present there is no mutual consensus as 

Figure 1.3: An example of a bent-core molecule. 
to the existence of the biaxial nematic phase, however the attention has been shifted to non-standardmesogens,suchasbent-coremoleculesandflexibledimers. Theproposition thatthe biaxial nematic phasemightbe field-stabilized [34] motivatesourwork presented in chapter 3. 
Bent-core mesogens are V-shaped molecules with a wide opening angle (typically 120 ∼ 140◦) (e.g. see Fig. 1.3), which received interest after their synthetization by Matsunaga and Matsuzaki in the 1990’s [35]. In these systems several smectic phases have been discovered [36] and also complex supramolecular structures [37], some discoveries point to existence of chiral smectic layers [38–40], which is peculiar for achiral molecules. In general, bent or bendable mesogens in appear to exhibit chiral symmetry breaking and ambidextrouschirality (spontaneousformation ofmultipledomainsofeither handedness) through conformation changes or forming of superstructures [41], such as mesoscopic helical and heliconical shapes, as first predicted by Dozov [42] and evidenced by Monte Carlo simulations [40, 43]. Other examples include spontaneous breaking of chiral and tetrahedratic symmetry in some ferrocene mesogens, which are allowed to bend due to the cyclopentadiene rings’ freedom to rotate with respect to each other [44]. Recently, in flexible dimers (chainstick-like mesogens with a flexible separator, which can attain bent shapes and also are achiral) there have been reports of an entirely new phase, dubbed the twist-bend nematic phase [45–49], in which the molecules form heliconical structures of either handedness. As evidenced by freeze-fracture experiments by the group from University of Colorado in Boulder [50], these structures have pitch on the nanometer scale, spanning only several molecules, unlike in other chiral liquid crystals, such as the cholesteric, where the pitch of the helix is hundreds to thousands of molecules long. In some sense these discoveries violate the usual intuition that chirality is necessarily intrinsic. We explore the phenomenon of spontaneous breaking of chirality in chapters 4 and 5. 
Symmetry considerations involving bent-core systems lead to a conclusion that a tetrahedratic order parameter should be included in their description [51, 52], which hints that phases of tetrahedratic symmetry may also be encountered. There are indications that this is indeed the case [53], furthermore it appears possible to induce tetrahedratic order by an external field in bent-core systems [54, 55]. The tetrahedratic phase T is optically isotropic (easily mistaken for the isotropic phase) and has the global symmetry of the point group Td, which is the symmetry group of the achiral regular tetrahedron (Fig. 1.4b). An achiral regular tetrahedron has four threefold axes, which pass through its vertices and the center, and three twofold axes, which lie on the midpoints of opposing edges. There are six mirror planes formed by two threefold axes and six mirror reflections combined with a rotation by 90◦. Together with the identity Td is of the order 24. The inclusion of tetrahedratic symmetry is discussed in chapter 4. 


1.3 Introduction to order parameters in liquid crystals 
The long history of active liquid crystal research has driven the development of a vast multitude of theoretical and computational approaches to the subject, much of which were applied from theories of phase transitions existing in statistical physics. While the interactions between mesogens are quite complex and calculating their approximations has grown into separate fields of research, much of this intrinsic complexity can be sacrificed when one isolates the crucial properties of the considered molecules. For instance, in uniaxial nematics it is sufficient to consider a unit vector cˆi, parametrized accordingly, to describe the orientation of the i-th molecule. (We remember that the translational degrees of freedom are free and at this point we do not need their explicit parametrization.) The normalized ensemble average
 cˆi 
=ˆn, (1.1)
| cˆ | 
is called the director and represents the average direction of long-range ordering of the molecules in the sample. An appropriate order parameter which takes into account how well, on average, the molecules are aligned with the director is defined as (see e.g. [2]): 
S =1  3cos 2θi − 1 , (1.2)
2
which is the ensemble average of the second Legendre polynomial P2(cosθi). One should notice that since the molecule (on average) looks the same when reflected along cˆi, the vector −cˆi is an equally good representation of the orientation of molecule i, thus a choice of cˆi · nn = cosθi is not a proper order parameter, as it averages to zero, due to lack of long range polar order. 
A thermotropic nematic above the clearing temperature forms an optically isotropic liquid, i.e. a phase of uniform spherical symmetry of the rotation group O(3). In such phase there is no preferred direction and therefore nˆis not a well-defined vector. The symmetry enables us to choose any frame of reference, and the average (1.2) yields zero. When the clearing point is approached with lowering temperature, the vector nˆis established by spontaneous symmetry breaking. If one knows nci,i =1, 2, ··· ,N, the instantaneous director and S can be computed. S becomes non-zero (in fact abruptly for uniaxial nematics) and continues to rise as the temperature is lowered. In real systems we eventually reach another phase transition, e.g. to a smectic phase or crystalline solid, dependingoncompound. Inanabstractsense, S would reachoneifthemolecules aligned precisely along the director. In experiment, S can be measured directly by e.g. NMR, birefringence or Raman scattering [2]. 
Calculating order parameters like (1.2) involves computing non-trivial integrals. For a general microscopic property A: 
  
A = Z−1 Tr A({Ωˆi,ηi})e−βH({Ωˆi,ri}), (1.3) 
{Ωˆi,ri}
  
−βH({Ωˆi,ri})
Z =Tr e , (1.4) 
{Ωˆi,ri}
where β =(kBT )−1 and H({Ωi,ri}) is the Hamiltonian of the system, which depends on the considered model (e.g. the Maier-Saupe Hamiltonian (1.8), discussed later). In the above, we have assumed that the kinetic part can be integrated separately. Thus, the trace is performed over the remaining rotational Ωˆi and translational ri degrees of freedom. In the case of rotational degrees of freedom, the trace is understood as SO(3)-invariant integration, therefore calculating (1.3) amounts to evaluating two 3Ndimensional integrals. Mean-field theory (section 1.7) and Markov chain Monte Carlo (section 1.8) provide methods for approximating (1.3), allowing one to escape the infea<BR>sibility of calculating such a large-dimensional integral. 

1.4 Dispersion forces 
The dispersion interaction, known in chemical literature as the London dispersion forces, first described by Fritz London (1900-1954) [56], is a type of Van der Waals interaction encountered in the second order of the quantum perturbation theory of interaction between two molecules at a large distance Rn12 [57]. The long-range character allows one to treat the interaction as a perturbation to the Hamiltonian of the free particle. The perturbative contribution to the energy of the pair ground state from dispersion interactions reads [58]: 
(1) (2) (1) (2)
 2
1 ψψ|H'|ψr ψs
00
udisp. = − , (1.5)
(1) (1) (2) (2)
2
Wr − W + Ws − W
=0,s=0
r��00 
(k)(k)(k)
where ψis the ground state of k-th molecule, while ψq are its excited states. W
00 and Wq (k) are energies of the corresponding state. The perturbation H' is simply the net electrostatic potential of interaction between all the charges residing at both molecules: 
(1) (2)
 ee
1 ij
H' = . (1.6)2rij
ij 
Subsequently, H' can be expanded in a multipoleseries and following premises depending on the considered molecules, some terms are discarded, others treated as dominating. Note that in the above expression both molecules are in excited states. Intuitively, the interactionsareduetothefactthatthermalfluctuations of thechargedistributiononone molecule result in a momentary excited multipolar state (dipole, quadrupole, octupole, etc.) which in turn induces an excited multipolar moment on the other molecule. Hence, the contribution is due to coupling between induced multipoles. 
In nematogens, which are mostly large, organic molecules, the dispersion term dominates overotherelectrodynamiceffectsbecauseitisproportionalto thevolume of themolecule. Other important contributions include electrostatic interactions due to static multipoles, such an electric dipole in bent-core molecules (see e.g. [36]). The dispersion interaction is easier to account for in calculations than purely electrostatic forces because of its effectively short range. In many cases the interaction can be constrained to the nearest neighbors without much sacrifice. However, it should be noted that the main cause for nematic ordering is excluded-volume (steric) effects, although the leading terms, in most cases, turn out to be similar to the leading dispersion terms. For example, dispersion and steric terms alike, in the case of the uniaxial, nematic yield interaction proportional to the second Legendre polynomial in the cosine of the intermolecular angle, as considered for dispersion interaction by Maier and Saupe [59–61] and as found by Onsager for the steric case [62]. Excluded volume effects produce phase transitions where the analogue of temperature is the inverse concentration 1 . To account for the temperature dependence 
ρ 
of the order parameters, thermotropic nematogens are studied mainly by dispersion-like models. For notable examples, dispersion interactions have been considered for uniaxial nematics by Maier and Saupe [59–61] (see also subsection 1.5.1), for biaxial nematics by Luckhurst et al. [63, 64], and for chiral symmetry breaking in bent-core molecules by Longa, Pająk and Wydro [1]. We will return to the last two models in chapters 3 and 4. 

1.5 Historical results 
1.5.1 Maier-Saupe model of thermotropic nematics (lattice case) 
An important result in the context of the presented work is the one of Maier and Saupe presented in a series of papers in 1958-1960 [59–61]. It is a specific application of the mean-field method, which we present in section 1.7, to a microscopic model of a uniaxial nematic, which allows for describing the thermotropic nematic phase transition. Maier and Saupe’s basic assumption was that the pair interaction potential between molecules i and j is due to point dispersion forces. The pair interaction potential is expanded in a multipole series in the the representation of spherical harmonics. By request, only uniaxially symmetric terms, i.e. with m =0, enter the expansion, which makes it sufficient to consider an expansion in Legendre polynomials. Therefore one obtains: 
α2
V (rij,θij)= − [P2(cosθij)+ ... ] , (1.7)
6
rij 
where rij is the relative distance between molecules i and j and θij is the relative angle between axes nci and ncj (as in Fig. 1.5). P1 and other odd terms are dropped since the interaction should be symmetric with respect to nci(j) →−nci(j), as polar ordering is being neglected. Inthegeneralcasethemean-fieldcalculationinvolvesanapproximationofthe two-point spatial correlation function by a step function, which is zero for rij ∈ [0,r0) and constant otherwise, to model the hard-core repulsion between molecules at close distances smaller than r0. However, if the molecules occupy sites of a simple cubic lattice and interact only with their nearest neighbors, one can put rij = a in (1.7), a being the lattice constant, thus lifting the need to integrate over positional coordinates. Neglecting higher-order terms, while taking α2/a6 ≡ α0, one arrives at: 

V (cosθij)= −α0P2(cosθij), 
(1.8)
1 
H = ij V (cosθij),
2 
which represents the Maier-Saupe l =2 minimal coupling model on the lattice. If we define the dimensionless temperature as t = β−1 = kBT/α0 (see e.g. (3.1)), the pair potential reads: 
1 
V (cosθij)= −P2(cosθij)= − (3cos 2θij − 1). (1.9)
2
Now assume that the system can be described by the single-particle distribution function f(θ) fortheorientationofmoleculeswithrespecttothedirector. Thisallowsonetowrite a single-particle average of the pair interaction potential, called the effective potential: 
Vef (cosθi)= dθjf(θj)V (cosθij). (1.10) 
Furthermore, f is expanded similarly as V in terms of Legendre polynomials: 
∞ 
f(cosθk)= f(nck · n
)= Pl Pl(nck · n
), (1.11) 
l=0 
where k numbers an arbitrary particle. Substituting (1.9) and (1.11) into (1.10): 
� ∞ 
Vef (cosθi)= Vef (nci · n
)= − dθj Pl P2(nci · ncj)Pl(nck · n
) (1.12) l=0 
and observing the orthogonality laws of Legendre polynomials, we arrive at: 
Vef (cosθi)= − Pl P2(nci · n
)= − Pl P2(cosθi). (1.13) 
In mean-field theory, (1.13) is used to construct the so-called equilibrium mean-field single-particle distribution, which can be used to calculate averages: 
Peq(cosθ)= Z−1 exp(−βVeff (cosθ)), (1.14) 
Z = dθexp(−βVeff (cosθ)), (1.15) 
P2 = dθP2(cosθ)Peq(cosθ). (1.16) 
Because P2 appears in (1.14) and (1.15), this set of equations is self consistent. We recognizethat P2 ispreciselytheparameter S asdefinedin(1.2). Forthesakeofbrevity we only indicate that these equations can be solved by numerical methods. The phase transition is of first order at Tc ≈ 0.22 with a sudden jump of S from 0 to S(Tc) ≈ 0.44. 

1.5.2 Lebwohl-Lasher lattice model 
The model of central importance as background to the presented work is the Lebwohl-Lasher lattice model of the uniaxial nematic. In 1972 Lebwohl and Lasher performed pioneering Monte Carlo simulations of a nematic in a lattice setup [65]. The setup consisted of molecules placed in the nodes of a 3-dimensional simple cubic lattice, with interactions restricted to the nearest neighbors and periodic boundary conditions, while the pair potential was adapted from (1.9): 
Vij = −EP2(cosθij), 1 (1.17)
H = Vij. 
2 
(i,j) 
The summation runs over nearest neighbors on the lattice, i.e. for each i summation over j runs over six neighboring particles. The authors considered 10 × 10 × 10 lattices far from the transition and 20×20×20 in the transition vicinity and considered sample sizes between 2000 and 8000 to discover a discontinuity of the mean energy, indicating a first-order phase transition at E/kBT =0.850 ±0.005 with a jump of the order parameter S of ΔS =0.33 ± 0.04. The results have been subsequently refined by Fabbrio and Zannoni for 30 × 30 × 30 lattices, who observed the phase transition at E/kBT =0.8903 ± 0.0005 [66]. Due to its simplicity and practical significance, Lebwohl and Lasher’s approach is often treated as a template in performing simulations of bulk liquid crystals. 


1.6 Alignment tensor 
In the general case, where apart from uniaxial ordering along director nˆ, the phase exhibits biaxial ordering (and two secondary biaxial directors {ˆl, mˆ} are established), it is concisely described by one tensorial parameter, called the alignment tensor: 
1 
Q = S(n 
⊗ n 
− I)+ T (nl − nm). (1.18)
l ⊗ nm ⊗ n
3 
Here S is the familiar uniaxial parameter, as defined in (1.2), while T is the biaxial order parameter. According to a definition by Straley [67]: 
T = sin2(θ)cos(2φ) , (1.19) 
where {φ, θ} are two of the three Euler angles, defined by consecutive rotations: Rz(φ)Ry(θ)Rz(ψ). One recognizes that the function being averaged on the right-hand side of (1.19) is (up to a normalization constant) the real part of the spherical harmonic Y l (θ, φ) for l =2,m =2, Y22(θ, φ). As we later discuss, in general the l =2
m
representation of SO(3) proves very useful in the context of biaxiality. 
The tensor (1.18) is traceless, symmetric and diagonal in the director basis {nl, nm, n
}, where it reads:
⎞
⎛ 
Q = 

⎜⎜⎝ 

−S 3 + T 00 0 −S 
3 
− T 

0 

⎟⎟⎠

. 

(1.20) 

2
00 S
3 
From the above form a few properties of Q can be inferred. For instance, if S �0
= and T =0, i.e. in the uniaxial nematic phase, Q is degenerate in the {nl, nm} plane, while the non-degenerate eigenvalue corresponds to the uniaxial director n
. Conversely, in the biaxial phase T �0 and all the eigenvalues are different. Because the tensor is
= traceless by definition, in the isotropic phase it equals to zero. Furthermore, the sign of the largest-modulus eigenvalue of Q informs us whether the uniaxiality is of the prolate (rodlike) (+) or oblate (discotic) (−) variety. 
It is often more practical to deal with Q directly than with S and T , as in phenomenological Oseen-Frank [68, 69] or Landau-de Gennes [2] theories. For instance, in the latter case, Q enters the expansion of the Landau free energy in the form of absolute rotational invariants Tr(Q2) and Tr(Q3). As it turns out, these are the only invariants without involving derivatives of Q, as from matrix algebra it follows that all invariants of higher order can be expressed by Tr(Q2) and Tr(Q3) [70]. These invariants hold a relation: Tr(Q3)2 ≤ 1 Tr(Q2)3 [70], which can be expressed with one condition w2 ≤ 1 [71, 72],
6 
where: 
√ Tr(Q3) 
w =6 . (1.21)
(Tr(Q2))3/2 For the uniaxial case w = ±1 and the sign distinguishes the prolate (+) and oblate (−) 
2
varieties. For the biaxial phase w< 1, while the case w =0 corresponds to the case of “maximum biaxiality”, when the eigenvalues of Q are {ω, 0, −ω} (in any order). 

1.7 Mean-field theory for lattice models 
As argued in section 1.3, an exact evaluation of the statistical average (1.3) is impractical in macroscopic systems. We present the mean-field method of producing an approximation for the equilibrium particle distribution, where the total number of degrees of freedom in the ensemble average is reduced by a factor of N. First, let us recall from statistical mechanics that the equilibrium state of the system in the canonical ensemble is given by minimum of the inequilibrium Helmholtz free energy functional [73, 74]: 
Fneq [P ]= U [P ] − TS [P ] , (1.22) 
where P = P (x1,x2,x3, ··· ,xN ) is the N-particle probability of finding the system in the state such that particle 1 is in the state x1, particle 2 is in the state x2, etc. We assume that the translational degrees of freedom are fixed and the N molecules reside on lattice sites, therefore the state xi of the i-th particle is parametrized by its orientation Ωˆi (xi ≡ Ωˆi). A trace over all possible states, Tr , is therefore understood as an SO(3)
{Ωˆi}
invariant integral over Ωˆi (e.g. for the parametrization of Ωˆi by Euler angles {αi,βi,γi}, such that the rotation matrix carrying from the laboratory reference frame to Ωˆi is R(Ωˆi)= Rzˆ(γ)Ryˆ(β)Rzˆ(α), where Ri(η) is a proper rotation about an angle η around 
 2π  π  2π
1
axis i, the trace reads Tr ≡ dαsinβdβdγ). At this point we assume that 
8π200 0 
{Ωˆi}
no continuous conformational degrees of freedom are present. If discrete conformational degrees of freedom are to be included, the trace needs to be supplemented by a sum over discrete conformational states. (An example is the chirality degree of freedom, introduced in section 4.2.) The equilibrium free energy is given by: 
Feq = Fneq [Peq] , (1.23) 
where Peq satisfies the conditional minimum equation: 
δ 
Fneq [P ] − λ Tr P (Ωˆ1 ··· ΩˆN ) − 1 =0. (1.24)
' 
δP {Ωˆi} 
We account for the constraint that P should be normalized through the Lagrange multiplier λ. (1.22) can be written explicitly using Shannon entropy: 
Fneq [P ]= Tr H(Ωˆ1 ··· ΩˆN )P (Ωˆ1 ··· ΩˆN )+ kBT Tr P (Ωˆ1 ··· ΩˆN )logP (Ωˆ1 ··· ΩˆN ) . 
{Ωˆi}{Ωˆi} 
(1.25) Integration over momenta is not taken into account, as for the models we are going to consider the kinetic part of the Hamiltonian gives an additive constant to the free energy. Generally, we assume that the Hamiltonian is a sum of pair interactions and single-particle potential energy (e.g. interaction with external field): 
NN
1 
H(Ωˆ1 ··· ΩˆN )= V (Ωˆi, Ωˆj) − V (Ωˆi). (1.26)
2 
(i,j) i 
Summation over i, j means that nearest-neighbor interaction is assumed, i.e. that each molecule i interacts only with its z neighbors. z is the coordination number, which depends on the lattice structure (e.g. z =6 for a simple cubic lattice with periodic boundary conditions). 
The main assumption of the mean-field approximation is that P (Ωˆ1, Ωˆ2, ··· , ΩˆN ) factorizes into identical single-particle distribution functions P (ˆ
Ω): 
P (Ωˆ1, Ωˆ2, ··· , ΩˆN )= P (Ωˆ1) · P (Ωˆ2) ·· ··· P (ΩˆN ). (1.27) 
With this assumption, (1.25) reduces to: 
N
1 
Fneq [P ]= Tr P (Ωˆi)V (Ωˆi, Ωˆj)P (Ωˆj)
2 {Ωˆi,Ωˆj }
(i,j) 
(1.28)
NN − Tr P (Ωˆi)V (Ωˆi)+ kBT Tr P (Ωˆi)logP (Ωˆi) . 
{Ωˆi}{Ωˆi}
ii 
Now solving the minimum condition (1.24) is a matter of applying the rules of functional calculus and algebra. The Lagrange multiplier λ is calculated from the normalization condition. The resulting equilibrium particle distribution function reads: 
−β(Veff (ˆΩ))
Ω)+V (ˆ
(ˆ
PeqΩ) = Z−1 e , (1.29)
eq 
Ω)+V (ˆ
Zeq =Tr e −β(Veff (ˆΩ)) , (1.30) {Ωˆ} 
where the “effective potential” is defined as: 
z 
ˆ
Veff (Ωˆj)= Tr Peq(Ωˆi)V (Ωˆi, Ωj) . (1.31) 
{Ωˆi}
i=j 
Furthermore, 
Feq = − N Tr Peq(ˆΩ) − NkBT logZeq. (1.32)
1 Ω)Veff (ˆ
2 {Ωˆ} 
The interpretation of (1.31) is that instead of involving the interactions of particle j with every of the other N − 1 particles explicitly, we treat j as if it were immersed in a field (Veff ) of interactions averaged over all the other particles, hence the term “mean-field” theory. The equations (1.29-1.31) are self-consistent, i.e. one generally cannot write a solution in closed form and they should be treated by e.g. iteration to convergence. 
The mean-field approach has the advantage of being relatively simple and computationally inexpensive. However, the transition temperatures are overestimated with respect to experiment and simulations, due to the assumption (1.27), i.e. that many-particle correlations are neglected. Because of this, entropy is underestimated, which enables the potential energy to dominate to higher temperatures than in reality. Another property of mean-field theory is that it produces Landau critical exponents regardless of the Hamiltonian in question (and therefore its applicability in studying critical behavior is limited). 

1.8 Monte Carlo simulations 
1.8.1 Metropolis algorithm 
To calculate the average (1.3), we also use the Markov chain Monte Carlo method, which in the modern version was described by Stanisław Ulam and Nicholas Metropolis in the 1940’s [75]. The name, owing to Ulam, is purportedly derived from the Monte Carlo casinos, in a reference of the method to gambling. The basics of the method stem from the Central Limit Theorem and sampling the equilibrium probability distribution by generating a Markov chain. 
The Central Limit Theorem, applied to the canonical ensemble, states that the ndimensional integral: 
−βH(x1,...,xn)
e
A = dx1 ··· dxnA(x1, ··· ,xn) (1.33)
Z 
can be replaced by an average over N samples {xi}: 
N→∞1 
A → A(x1, ··· ,xn), (1.34)
N 
{xi}:p({xi}) 
where the samples {xi} are drawn from: 
−βH(x1,...,xn)
p({xi})= Z−1 e . (1.35) 
An algorithm which allows sampling of p({xi}) is given by Metropolis [76]. The prescription amounts to designing an ergodic Markov process that asymptotically reaches p({xi}) as the stationary distribution. This is done by observing that such a process needs to conform to detailed balance: 
p({xi})σ({xi}→{xj})= σ({xj})p({xj}→{xi}), (1.36) 
where σ({xi}→{xj}) is the probability of transition of the system from state i to j. Detailed balance secures that the transition is reversible, i.e. that the considered system is in equilibrium. The equilibrium probability distribution is given by (1.35), so rewriting (1.36): 
p({xi}→{xj})e−βH({xi}) 
−βΔEij
= =e , (1.37)
p({xj}→{xi})e−βH({xj }) 
where ΔEij is the difference of energy between the two equilibrium states. The prescription due to Metropolis is that we can achieve the relative probability (1.37) in a simulation of the system by following an algorithm: 
1. 
Starting with the system in the state i, make a small random change in i, thus creating the “trial” state j and calculate the energy difference ΔEij. 

2. 
If ΔEij ≤ 0, accept j as the new state of the system unconditionally. 

3. 
If ΔEij > 0, accept j as the new state of the system with the probability e−βΔEij . 

4. 
Repeat 1-3 for a suitable sample size. 


This algorithm produces a random walk in the state space of the system. For a good convergence to the real distribution p({xi}), the space of equilibrium states needs to be explored quite extensively. The question is how large the change to i in pt. 1, i.e. the random-walk step should be. Very large changes will almost never be accepted, while very small ones will have the system stuck in one area of the state space for a very long time. Generally, a rule of thumb is assumed that the step should be of such magnitude, that the fraction of accepted moves is around 0.4 ÷ 0.5. In most cases the changes will be small enough so that states i and j will be substantially correlated, thus making the random walk non-Markovian. To obtain a good Markov chain, one needs to discard a number of intermediate states until no correlation is measured. The final result is a Monte Carlo sample of the statistical ensemble, i.e. a set of replicas of the system in equilibrium. On the other hand, ergodicity permits us to regard the replicas as snapshots of a single system undergoing an equilibrium process. 
There remains the problem of the choice of the initial state i, since an arbitrary state will likely be out of equilibrium. While the Metropolis algorithm is designed to produce an approximation for the equilibrium distribution function, it can be used to find an equilibrium process through random walk in state space, which drifts towards an energy minimum, in a non-equilibrium process called thermalization. Once thermalization is complete, sampling of equilibrium states can begin. A question open to the scientist is to make sure that the system is not stuck in a metastable state. 

1.8.2 Rotational degrees of freedom 
In the presented work the simulations are performed in the spirit of Lebwohl and Lasher (see section 1.5.2), i.e. the Monte Carlo degrees of freedom are limited to the orientation of molecules (in later chapters we also consider chirality). There are several parametrizations of rotations in 3D space, e.g. Euler angles, which describe the subsequent rotations of the frame of reference around axes zˆ,xˆand zˆ. The computational disadvantage of this and many other methods is the need of evaluating trigonometric functions. In Monte Carlo simulations such evaluations are performed for each particle in each attempted move, which can turn out to be substantially expensive. For this reason we parametrize rotations using quaternions, which reduces the problem to evaluating second-order polynomials [77]. 
The parametrization follows from the fact that the group SU(2) of unitary 2×2 matrices of determinant 1: 
u1 u2 
U = , (1.38)
∗∗
−uu
21 
U† = U−1 , (1.39) 
22
detU = |u1| + |u2| =1, (1.40) 
isadoubleuniversalcoverof SO(3), thegroupof 3×3 orthogonalmatricesofdeterminant 
1. On the other hand, SU(2) is isomorphic with the group of quaternions with norm 1: 
q0 + iq3 iq1 − q2 
U = , (1.41)iq1 + q2 q0 − iq3 
q = q0 + q1i + q2j + q3k, (1.42) 
2222
|q|2 = detU = q0 + q1 + q2 + q3 =1, (1.43) 
where q is a quaternion. Thus, the space of rotations in 3D space is replaced by a set of points on the unit 4D sphere. The mapping is given by expressing the rotation matrix 
in terms of q (see e.g. [77]): <BR><BR>
⎞
⎛ 
2222
q0 + q1 − q2 − q2(q1q2 − q0q3) 2(q1q3 + q0q2)
3 
R = 

⎜⎜⎝

2(q2q1 + q0q3) 
q

2222 
0 − q1 + q2 − q
2(q2q3 − q0q1) 
⎟⎟⎠

. 

(1.44)

3 
2222
2(q3q1 − q0q2) 2(q3q2 + q0q1) q− q− q2 + q
01 3 
Note that R does not change if q →−q, which reflects the fact of double covering. The matrix (1.44) gives us the reference frame rotated according to q. As stated before, when comparing the rotation matrix for Euler angles and (1.44), it is readily seen that the quaternion version is computationally cheaper to evaluate. 
To make use of quaternions in Monte Carlo, we need a means of generating random rotations expressed in quaternions. The parametrization (1.41-1.44) simplifies this task by reducing it to drawing points from an uniform distribution on a 4-sphere, because the invariant measure for SU(2)-invariant integration is uniform for quaternions, i.e. 
22
Tr ≡ δ(q0 + ··· + q3 − 1)dq0 . . . dq3 for the canonical ensemble average (1.3). The 
{Ωˆi}
algorithm used in the presented work is due to Marsaglia [78]. First, we draw two independent points from the unit circle: 
p1 =(y1,y2),p2 =(y3,y4),yi ∈ [−1, 1], (1.45) 
2 22 
r = y1 + y2 ≤ 1, (1.46)
1 
2 22 
r = y3 + y4 ≤ 1. (1.47)
2 
If so, the point: 

22
(y1,y2,y3 (1 − r)/r22,y4 (1 − r)/r2) (1.48)
112
is uniformly distributed on the unit 4-sphere. Thus, we simply need to draw a quaternion with accordance to (1.48) to produce a random rotation in 3D space. 
However, if one makes a trial move on a molecule by simply replacing the quaternion describing its rotation by a new, random one, the relative change will be most likely quite large and many moves will be rejected. Therefore, we perform the trial move by slightly altering the present orientation of the molecule. Considering that q is the quaternion for the molecule in question, we draw a random quaternion Δq according to (1.48). Then, 
'
the new quaternion q is: 
q + rΔq
' 
q = , (1.49)
|q + rΔq|where 0 <r ≤ 1. The parameter r can be adjusted to ensure proper acceptance ratio of trial Monte Carlo moves. To calculate the difference in energy, one can rewrite the Hamiltonian using quaternions and do away with (1.44) completely during the calcula<BR>tion. We, however, found that (1.44) bears an infinitesimal cost worth the advantage of retaining the original expression in terms of e.g. Cartesian tensors built on the molecular basis {eˆx,eˆy,eˆz} given by R. 


1.8.3 Parallel sampling and pseudorandom number generation 
Once we have an equilibrium state at our disposal, we can take advantage of parallel computing to either increase the sample size or reduce the time needed for the simulation tocomplete. Thisisperformedformallybybranchingtherandomwalkinstatespaceinto k random walks, which can explore the state space independently. If the random walks are non-correlated, the generated sets of equilibrium states can be treated as one sample of the statistical ensemble. Computationally, this is performed by first thermalizing the system into equilibrium and then taking k replicas, which are considered as initial states of k independent Monte Carlo simulations in separate threads. The results are aggregated at the end. 
The principal problem of this method is ensuring that the resulting Markov chains are non-correlated, which amounts to securing that each of the threads has its own, unique sequence of pseudorandom numbers. This issue is studied extensively and has been addressed recently for the Mersenne Twister algorithm on graphic processors by Podlozhnyuk [79] and earlier for linear congruential generators by Durst [80]. A method involving cellular automata has been proposed by Hortensius et al. in 1989 [81]. In the presented work we use a relatively new pseudorandom number generator family invented by Marsaglia in 1994 [82, 83] popularly referred to as “multiply-with-carry”. The sequence of the generator is defined as: 
xn+1 =(axn + cn) mod b (1.50) 
where a is a multiplier and b is a modulus, which is typically 232 for 32-bit and 264 for 64-bit processor architectures. cn is the carry, defined as: 
axn−1 + cn−1 
cn = lJ (1.51)
b 
The multiplier a needs to be chosen such that ab − 1 is a safe prime, i.e. both ab − 1 and (ab − 1)/2 − 1 are prime. The period of (1.50) is of the order of b, thus greatly improved by utilizing the 64-bit architecture. The advantage of (1.50) is that one can produce infinitely many generators for a given b. Thus, at the expense of a small period, when compared to e.g. that of Mersenne Twister, equal to 219937 − 1, we can generate non-correlated sequences of pseudorandom numbers by specifying different multipliers and seeds to separate threads, instead of sharing one generator, which comes with implementation issues. In our software, we use the multipliers and seeds for modulus 232 generated by Gratton [84]. 


1.9 Purpose and plan of thesis 
The aim of the present thesis is to study some of the remarkable effects observed in systems such as bent-core and ferrocene mesogens, as well as flexible dimers, including field-stabilized biaxiality, homogeneous chirality and twisted states of homogeneous and ambidextrous chirality, through a lattice dispersion model with the help of the methods cited in this chapter – primarily by Monte Carlo simulation and, to a lesser extent, mean-field theory. The dispersion model being studied is an extension on models previously considered in literature, of which all are generalizations of the Lebwohl-Lasher lattice model, discussed in section 1.5.2. The definition of this general model, along with descriptions of the respective terms and their purposes is given in Chapter 2, after a brief introduction to the formalism of irreducible tensors and symmetrization. Since this model is quite broad in scope and embraces six independent parameters, our aim is constrained to investigating several limit cases of the general model. 
In Chapter 3 the case of a dispersion model of biaxial nematics in an external field is considered. The interaction potential (3.14) is based upon the model by Luckhurst and Romano [64], augmented by a term which is quadratic in field (see (3.12)). The field considered is a generalized vector field, thus the effects apply (almost) equally to the magnetic and electric cases. The field-interaction potential promotes alignment of either the molecular axis of largest anisotropy or the largest molecular face, depending on the sign of the anisotropy of the generalized (dielectric or diamagnetic) molecular polarizability. Furthermore,theLuckhurst-Romanomodel,expressedin(3.1),introduces 
 
a parameter for controlling the biaxiality of molecules, λ, which ranges from 0 to3/2. 
The values correspond to prolate (rodlike) and oblate (discotic) shapes for λ< 

1√ 
and λ> 

6 
1√
respectively, while λ = 

6 
1√
is a highly-symmetric value, for which the model 

exhibits a direct phase transition from the isotropic phase to a biaxial phase of maximum 
biaxiality. Therefore, we study the model with external field considering different values 

of λ (λ =0.3,λ =0.5 and λ 

= 

6 
1√
) and the sign of the anisotropy of the molecular 

polarizability (positive and negative), while observing the effects of varying magnitude of the external field. 
The liquid crystalline phases originally known to exist for the Luckhurst-Romano model, presented in the phase diagram in Fig. 3.3, are the isotropic (I), uniaxial nematic prolate and oblate (NU+ and NU−) and biaxial nematic (NB). After an introduction of order parameters and the details of the theoretical methods used (section 3.4), as a test of out methods we present results for this model in section 3.5. In section 3.6 we proceed to presentation of results for the model with external field, starting with an analysis of the critical and tricritical points of the isotropic-uniaxial nematic transition in section 3.6.1. We then move on to the Monte Carlo results, compared with mean-field predictions, for λ =0.3, for which value the mean-field phase diagram in the (t, ΔEh2)-plane (t being the dimensionless temperature and h being the magnitude of the field) is compared with the energy fluctuation landscape obtained from simulation, in section 3.6.2. For λ =0.5 the mean-field phase diagram is obtained, presented in section 3.6.3. For the highly
symmetric value of λ = 

1√
the mean-field phase diagram is obtained and supplemented 
6 
by temperature scans of the order parameters, fluctuation of energy and susceptibilities of order parameters in section 3.6.4. 
In Chapter 4 the case of the biaxial nematic model with added tetrahedratic coupling (withoutthepresenceofafield)isinvestigated. Theinteractionpotential,firstpostulated in [1], consists of the quadrupolar coupling term, as in the Luckhurst-Romano model, and an octupolar coupling term, scaled by the coupling constant κ. The relationship between the D2h symmetry of the molecular quadrupolar moment and the Td symmetry ofthemolecularoctupolarmomentproducestwopossiblechiralmolecularconfigurations, distinguished by the parity of the molecular basis, which is considered as a binary degree 
28
of freedom. We study this model for τ =1 and a highly-symmetric value τ = with
15 
Monte Carlo simulations. For both cases there exist phases of tetrahedratic and chiral symmetry, while for the latter case a direct transition from the isotropic to the chiral phase exists. The results are presented in sections 4.4 and 4.5 respectively. 
The results contained in Chapter 5 are meant to be approached as preliminary. In this chapter we consider the model studied in Chapter 4 with an added term, which accounts for coupling of the molecular octupolar moment to the intermolecular lattice vectors. Asarguedinsection5.2,thiskindofcouplingproducestwistsbetweenmolecules on neighboring sites, as is promptly demonstrated in simulation on a one-dimensional chain. Subsequently, a two-dimensional case is considered in section 5.3, for which the chiral-symmetry-breaking phase transition is studied and two types of structures in the chiral phase are identified. The three-dimensional case is discussed in section 5.4 and as of the moment of completion of this thesis the results presented are still largely open to interpretation. Two types of structures are found and presented in Fig. 5.6. Words of conclusion and comment on further studies are found in the closing Chapter 6. 
Chapter 2 


Generalized dispersion model for 
bent-core (and related) systems 

Before we define the model studied here itself, we first provide an introduction to the tool which is central to the understanding and formulation of the former – the formalism of symmetrized irreducible spherical tensors. An explicit definition of the model in question is given subsequently in section 2.2. 
2.1 Irreducible tensors and symmetrization 
From the argument in section 1.4 it follows that dispersion interactions amount to cou<BR>pling of multipolar moments, which can be expressed in terms of Cartesian tensors. For dispersion models which possess global rotational symmetry this representation is however often unoptimal. A more practical approach is to express the tensors involved in the dispersion interaction in terms of spherical tensors, which obey much simpler transformation laws when undergoing rotations, i.e. they are always confined to the same fixed-l subspace of angular momentum [85]. 
Irreducible spherical tensors of rank n (the number of Cartesian indices), for angular momentum quantum number l (l =0, 1, ··· ,n − 1,n) and projection of angular momen
(l),n
tum quantum number m (m = −l, −l +1, ··· ,l − 1,l) are denoted as eˆm . They can be expressed within an ordinary Cartesian basis {eˆx,eˆy,eˆz}, which can be later on identified with the molecular basis or the director tripod, depending on our needs. Starting from tensors of rank 1: 
(1) (1) 1 
eˆ=ˆez,eˆ= �√ (ˆex ± ieˆy), (2.1)
0 ±1 
2
23 

one can recursively construct tensors of arbitrary rank with the help of Clebsch-Gordan coefficients: 
l − 11 
l
(l,δ),n 
(l−1),n−1 (1),1 
eˆ= 
eˆ⊗ eˆ. (2.2)
m 
m1 m2 
m
m1 m2
m1,m2 
Here δ is the “seniority index”, used to distinguish between linearly independent same-l tensors within the same Cartesian rank n, which appear for n ≥ 3. From now on we decide to drop the seniority index, except where necessary. Under a rotation R, the irreducible spherical tensor transforms through the Wigner rotation matrix: 
(l),n (l),n (l)
Reˆm = eˆ' D', (2.3)
mmm
'
m
��∗ �� 
(l)(l),n (l),n
D'= eˆ·Reˆ' , (2.4)
m
mm m
where the scalar product “·” is understood as a full contraction over all Cartesian indices. As one notices from the above expression, a rotation always leaves the tensor in the same l-subspace of angular momentum, thus mixing only tensors with different m. From the many remaining useful properties of irreducible spherical tensors, for the sake of brevity we only name orthonormality between tensors in the same-l subspace. For a comprehensive guide on the properties and construction of spherical tensors and Wigner matrices, cf. e.g. [86]. 
In the context of liquid crystal researchthe mainadvantage from irreducible tensors flows from symmetry-adaptation. Considering different subgroups of SO(3) (or, in general, O(3)), one can construct irreducible tensors which bear the symmetry of the individual subgroups. The obtained set of tensors gives us complete control over the symmetry in designing the interaction potential and also serves as a neat, symmetry-adapted basis for the alignment tensor. The process of symmetrization is straightforward. If the group in question is G, the G-symmetrized tensor is obtained using the prescription: 
1 
T(l),n l (l),n
= c O(g)ˆe. (2.5)
mm m
|G| 
g∈G 
Here |G| stands for the number of transformations in G and O(g) is the linear operator 
l
associatedwithtransformation g. Thecoefficients carechosentosecureorthonormality 
m 
of tensors with the same m. In the case of continuous groups, the sum needs to be replaced with SO(3)-invariant integration and |G| is the unit integration volume of the group. 
Having at hand the tensors (2.5), we easily produce the complementary symmetrized Wigner matrices, which we denote Δ(l) ' , following the notation of Mulder [87]: 
mm 
(l) = T(l),n (l),n
Δ(Ωˆ1) · T(Ωˆ2), (2.6)
mm ' mm ' 
where Ωˆ1 = {eˆx,eˆy,eˆz} is often identified with the director basis, and Ωˆ2 is some other frame of reference, commonly identified with the molecular basis. 
Symmetrized tensors form a complete basis for the l-subspace for a given symmetry group G (a set involving all G-symmetrized tensors for angular momentum l), therefore we have an expansion of the (G-symmetrized) identity operator: 
(l,δ) T(l,δ) Ω) ⊗ T(l,δ)
(ˆ(ˆ
I= Ω). (2.7)
G mm m,δ 
2.1.1 Symmetrized irreducible tensors for selected symmetry groups 
For the use in this chapter and for the remainder of this thesis we consider the leading symmetrized tensors for l =2 the D2h-symmetrized second rank tensors and for l =3 the Td-symmetrized third rank tensor. 
The non-vanishing l =2 symmetrized tensors for the biaxial symmetry group D2h, with symmetry axes coinciding with {eˆx,eˆy,eˆz}, are: 
   
T(2),2 0  =  3 2(ˆez  ⊗ ˆez  − 1 3 I),  (2.8)  

1
(2),2
T= √ (ˆex ⊗ eˆx − eˆy ⊗ eˆy) (2.9)
2 
2
√
(0),2
and the trivial isometric tensor eˆ=1/ 3I, where I =ˆex ⊗eˆx +ˆey ⊗eˆy +ˆez ⊗eˆz is the 
0 
Cartesian identity operator. One recognizes that, when Ωˆ= {ˆl, ˆn}, (2.8) and (2.9) 
m, ˆappear, up to a normalization constant, as components of the alignment tensor (1.18). Furthermore, the l =2 D2h-symmetrized Wigner matrix elements Δ(2) are, up to a 
m,m ' 
(2)
normalization constant, the functions identified by Straley [67], in particular S ∝ Δ0,0 
(2)
and T ∝ Δ.
0,2 
For the tetrahedratic group Td there are no non-vanishing symmetrized tensors of l< 3, apart from trivial isotropic tensors. Distinction by seniority of rank 3 tensors is only needed for tensors of angular momentum l =2, therefore it can be dropped altogether. 
(3),3
For n =3,l =3 there are 7 independent tensors eˆm , m = −3 ... 3, and depending on the choice of the symmetry axes, a symmetrized tensor will be a linear combination of l =3 tensors with different m. With the choice (that we consider later) that the twofold axes of symmetry coincide with {eˆx,eˆy,eˆz}, the only non-vanishing Td-symmetrized irreducible tensor is [1]: 
1
(3),3
T2  = √ 6 (ˆei,ˆej ,ˆek)∈π(ˆex,ˆey,ˆez) ˆei ⊗ ˆej ⊗ ˆek,  (2.10)  

where the sum runs over all permutations of {eˆx,eˆy,eˆz}. 


2.2 Definition of the model 
The most general model studied here is a lattice dispersion model, defined by the Hamiltonian: 
NN
1 
ˆˆˆˆ
H = VQ(Ωˆi, Ωj)+ VT (Ωˆi, Ωj)+ Vc(Ωˆi, Ωj)+ V T (Ωˆi, Ωj)+ ΔV (Ωˆi). (2.11)
c
2 
(i,j) i 
The sum runs over nearest neighbors, i.e. for each particle i the sum runs over z =6 neighboring particles j, while the N particles occupy sites of a simple cubic lattice with periodic boundary conditions (although two other options are considered in Chapter 5). Thearguments Ωˆi andΩˆj denoteorientationsoftheparticlesatlatticesitesi andj,which can be expressed equivalently by specifying three orthonormal vectors – the molecular 
ˆˆ
frame of reference (also: “molecular tripod”, “molecular basis”): Ωi = {aˆi,bi,cˆi} and ˆˆ
Ωj = {aˆj,bj,cˆj}. Note that the molecular frame of reference can be either righthanded or lefthanded, which will become important when we take a closer look at VT . The interaction terms are explained as follows: 
• The quadrupolar term: 
VQ(Ωˆi, Ωˆj)= −EQ(Ωˆi) · Q(Ωˆj) (2.12) 
couples molecular quadrupolar moments at sites i and j. The scalar product ’·’ is understood as a full contraction over Cartesian indices. The quadrupolar moment Q(Ωˆi) is a second-rank traceless tensor, which has uniaxial and biaxial parts, defined in terms of irreducible spherical tensors symmetrized with respect to the uniaxial (D∞h) and biaxial (D2h) symmetry groups (2.8-2.9), expressed in the molecular basis: 
√
(2) (2)
Q(Ωˆn)= T(Ωˆn)+ 2λT(Ωˆn),n = i, j. (2.13)
02 
The quadrupolar term VQ is responsible for uniaxial and biaxial nematic ordering. Its history, motivations and previous results are discussed in Chapter 3. 
• The octupolar term: 
(3) (3)
ˆ
VT (Ωˆi, Ωj)= −τET(Ωˆi) · T(Ωˆj) (2.14)
22 
couples molecular octupolar moments at sites i and j. T(3)(Ωˆi) is a third-rank 
2 
irreducible spherical tensor, symmetrized with respect to the tetrahedratic group Td, of the form given in (2.10), expressed in the molecular basis. Together with VQ, this term is shown to be responsible for emergence of phases of tetrahedratic symmetry and spontaneous breaking of chiral symmetry, recently observed experimentally. Chirality enters the model through parity (handedness) of the molecular basis: 
pi =ˆai · (ˆbi × cˆi)= ±1, (2.15) 
(3)(ˆ
because a parity switch is a symmetry of Q(Ωˆi), but not of T2 Ω). In this way, conformations of different chirality are accounted for as microscopic states parametrized by the parity degree of freedom. Further explanation is given in Chapter 4. 
• The intermolecular vector coupling terms: 
ˆ
Vc(Ωˆi, Ωj)= κE [ˆbαi · (ˆbβi × ˆbγi)Qαν(Ωˆi)Qβν (Ωˆj) 
(2.16) 
− ˆbαj · (ˆbβj × ˆbγj)Qαν(Ωˆj)Qβν(Ωˆi) (ˆrij)γ, 
(3) (3)
V T (Ωˆi, Ωˆj)= κT E [ˆbαi · (ˆbβi × ˆbγi)T (Ωˆi)T (Ωˆj)
c 2 ανµ2 βνµ
(2.17)
(3) (3)
− ˆbαj · (ˆbβj × ˆbγj)T (Ωˆj)T (Ωˆi) (ˆrij)γ ,
2 ανµ2 βνµ
Q(ˆ(3) (ˆ(3)(ˆˆ
whereα, β, γ, ν = {x, y, z}, Qρδ(ˆΩ) ,T Ω) = TΩ) ,bxi =
Ω) = 
2 αρδ2
ρδ αρδ 
aˆi, ˆbyi =ˆbi, ˆbzi = cˆi and Einstein summation convention is implied, couple the molecular quadrupolar (Vc) and octupolar (V T ) moments to unit vectors rˆij con-
c 
necting sites i and j. These terms should be considered because the antisymmetric tensor ˆbαn · (ˆbβn × ˆbγn),n = i, j can be constructed out of the other tensors: 
√ 
(2) (2) (3)
ˆ
bαn · (ˆbβn × ˆbγn)=2 2 T0,xµ T2,yνT2,µνz, (2.18) (x,y,z)∈π{α,β,γ} 
(T )
where the summation runs over cyclic permutations of {α, β, γ}. Terms Vc will always be generated in the presence of chirality and will lead to an inhomogeneous ground state of the model. (It should be noted that we imply the resulting inhomogenities to be of large characteristic lengths when compared to the lattice spacing.) Further explanation is given in Chapter 5. 
• 	The field-coupling term: 
ΔV (Ωˆi)= −EΔEhαhβQ(Ωˆi)αβ 	(2.19) 
couples the molecular quadrupolar moment Q(Ωˆi) to an external, general vector field nh. α, β = {x, y, z} and Einstein’s summation convention is implied in the above definition. EΔE is the anisotropy of the molecular (dielectric or diamagnetic) field polarizability, where ΔE = ±1, thus the anisotropy is taken equal to the coupling constant in the quadrupolar term VQ above or as a renormalizing factor to the squared field magnitude h2. The origin of this term is further explained in Chapter 3. 
The model (2.11), is designed to reproduce the most remarkable recent findings made in bent-core systems, ferrocene mesogens and flexible dimers, i.e. biaxial phases [27, 29, 30, 32, 88], spontaneous breaking of chiral symmetry, formation of modulated chiral struc<BR>tures [38–40, 43] including the twist-bend nematic phase [45–50], and tetrahedratic order [51, 52, 54, 55]. In this most general form there are seven independent parameters: the energy scaling parameter E, the biaxiality parameter λ, the tetrahedratic coupling constant τ , the intermolecular vector coupling constants κ and κT , the sign of the anisotropy of the molecular field polarizability ΔE and the magnitude of the external field h = |nh|. E can be eliminated from the list of free parameters by studying the model with respect to the dimensionless temperature t = kBT/E (which is equivalent to putting E =1). Furthermore, only one of κ and κT needs to be considered, as they generate identical ground states, therefore we choose to disregard κτ . Nevertheless, studying the full spectrum of the remaining five parameters is a large task in itself, but, more importantly, we might have trouble interpreting the results without first studying limit cases. Therefore, this thesis presents a study on several special cases of the general model (2.11). Namely: 
• 	In Chapter 3 the case of E> 0,λ > 0,h > 0,τ = κ =0, i.e. the l =2 model of biaxial nematics with external field coupling is studied in detail with Monte Carlo simulations and mean-field methods. This limit case, as we will show, produces effects such as field-induced biaxiality in the uniaxial nematic phase, as well as critical and tricritical points of the isotropic-uniaxial nematic transition as the field 
1
is increased. For a highly symmetric value of λ = √ an external field produces 
6 
an unusual effect of reducing the transition temperature and decreasing the level of biaxiality in a highly-symmetric biaxial phase. 
• 	In Chapter 4 the case of E> 0,λ > 0,τ > 0,h = κ =0, i.e. the model proposed in 
[1] is studied for selected values of λ and τ with Monte Carlo simulations. This case extends on the studies in [1] by giving full characteristics of homogeneously chiral 
28
and tetrahedratic structures. In addition, a highly-symmetric value of τ = is
15 
1
studied, for which at λ = √ there exists a direct transition from the isotropic to 
6 
the chiral phase. 
• 	In Chapter 5 selected Monte Carlo simulation results for the case of E> 0,λ > 0,τ > 0,κ > 0,h =0 are presented. This variant extends on the previous one by accounting for intermolecular vector coupling, which leads to spatially inhomogeneous phases. We present only preliminary result for this case for 2D and 3D systems. In both cases two distinct structures are identified in the chiral phase, described in more detail in sections 5.3 and 5.4. 
Chapter 3 


Dispersion model of biaxial 
nematics in an external field 

3.1 Introduction 
The possibility of stabilizing a thermotropic biaxial nematic phase in bent-core mesogens [27, 29] has invigorated the interest in this elusive phase, particularly in the possibility of magnetic or electric field-induced biaxiality [34]. In fact, there are several accounts of field-induced biaxiality in bent-core systems [89–92], apart from numerous results where magnetic or electric field effects produce biaxiality in other systems [93–95]. Although several theoretical descriptions exist [70, 96–98], along with various simulational studies [99–102], a thorough analysis of field-induced effects by a microscopic model appears to be missing. We will now consider a lattice dispersion model of a biaxial nematic liquid crystal interacting with an external vector field, which will be studied with mean-field andMonteCarlomethods. Thischaptercompletesthepreliminaryresultsfirstpublished in [103]. 

3.2 Dispersion models of biaxial nematics 
The idea of a biaxial nematic phase first occurred to Freiser in 1970 [21]. He asserted that for correct description of biaxial systems the l =2 subspace needs to be considered in the multipole expansion of pair dispersion interaction energy and predicted that a sequence of a first-order transition to the uniaxial nematic and second-order transition to the biaxial nematic should occur with lowering temperature. 
31 

In 1974 Straley [67] considered a generalization of the Maier-Saupe pair potential (1.9) constructed on D2h-symmetry-adapted functions for l =2. Assuming a system of bricklike objects of length L, breadth B and width W , using mean-field methods he produced the phase diagram for objects of varying breadth. It appears so that there exists a special 
√ 
value of B = LW , corresponding to maximum biaxiality, where the transition from the isotropic phase is directly into the biaxial phase. Furthermore, Straley noticed that the 
√√ 
models for rodlike (B< LW ) and platelike objects (B> LW ) can be mapped onto 
√ 
each other through a duality transformation. At B = LW the transformation carries 
√ 
the model into itself, therefore B = LW has been named the “self-dual” point. 
Similarly to Freiser, Luckhurst and Zannoni considered a l =2 expansion of the point dispersion potential in Wigner rotation matrices D(2) ' and developed a molecular field 
mm 
model [63]. Monte Carlo simulation for a l =2 model has been first performed by Luckhurst and Romano [64]. Their version introduced the parameter λ, which controls 
3
the biaxiality of molecules, closely related to Straley’s B. Namely, 0 ≥ λ ≤ , while 
2 
the boundary cases correspond to Maier-Saupe-like uniaxial models for long rods (λ =0) 
31
and platelets (λ = ). The case of maximum biaxiality is represented by λ = √ . The 
2 6 
phase diagram for varying λ has been obtained using mean-field theory and supported by Monte Carlo simulations by Biscarini et al. in 1995 [104]. Thorough Monte Carlo investigation of the self-dual point has been performed by Chiccoli [105]. 
The l =2 model family has been studied and reformulated extensively. Mulder explored 
(2)
the expansion in terms of symmetry-adapted functions Δmm ' with bifurcation theory [87]. A version of the model has been studied by Sonnet et al. [106] and Bisi et al. [107] for which an universal mean-field phase diagram has been obtained. Density functional and bifurcation theory has been analyzed by Longa et al. [108]. In parallel, comprehen<BR>sive Landau-de Gennes theory has been developed by Palffy-Muhoray and Dunmur [96], Allender [109], Gramsbergen et al. [70] and more recently by Allender and Longa [72]. Longa and Pająk further analyzed the link between phenomenological theory and the l =2 Luckhurst-Romano model [71]. 
3.2.1 Dispersion model of interacting quadrupoles 
The base for our model is the Luckhurst-Romano model expressed in the parametrization ofirreduciblesphericaltensors. Inthismodelitisassumedthatthedispersioninteraction is due to coupling of induced quadrupolar moments. In the general model (2.11), this model corresponds to the case where all parameters aside from λ and E are zero. The pair potential and Hamiltonian read: 
ˆ
VQ(Ωˆi, Ωj)= − EQ(Ωˆi) · Q(Ωˆj), 
N
1 (3.1)
ˆ
H = VQ(Ωˆi, Ωj),
2 
(i,j) 
ˆˆ
where E> 0 and Ωˆi = {aˆi,bi,cˆi} and Ωˆj = {aˆj,bj,cˆj} specify orientations of the molecular frames of reference with respect to the laboratory frame of reference {ˆy, ˆ
x, ˆz} for molecules i and j respectively (see Fig. 3.1a). The sum runs over nearest neighbors, i.e. for every i we take z =6 neighboring particles on a simple cubic lattice with periodic boundary conditions. The scalar product “·” is understood as a full contraction over all Cartesian indices. The quadrupolar molecular moment Q(Ωˆi) is a second-rank tensor, defined as: 
√
(2) (2)
Q(Ωˆi)= T(Ωˆi)+ λ 2T(Ωˆi). (3.2)
02 
3 (2)
λ ∈ [0, ] is a parameter which controls molecular biaxiality by combining T(Ωˆi)
20 and T(2)(Ωˆi) – the second-rank D∞h and D2h-symmetrized irreducible spherical tensors 
2 
for l =2, previously defined in (2.8-2.9), which here are constructed in the molecular basis Ωˆi: 
 
(2) 31 
T(Ωˆi)=(ˆci ⊗ cˆi − I),
0 
23 (3.3)
1 
T(2)2 (Ωˆi)= √ (ˆai ⊗ aˆi − ˆbi ⊗ ˆbi). 
2
As the interaction (3.1) is regarded as the leading part of the intermolecular dispersion interaction, Q(ˆ
Ω) can be interpreted as being proportional to the anisotropic part of the molecular (dielectric or diamagnetic) polarizability tensor α. In this case, the parameters E and λ are related to the diagonal elements of the diagonalized tensor α {αx,αy,αz}: 
)23 αx−αy
E = E0(2αz − αx − αy, λ = .
22αz−αx−αy 
If one expresses (3.1) using the symmetry-adapted Wigner functions (2.6) for l =2, by ˆˆ
substituting Ωˆ1 =Ωi and Ωˆ2 =Ωj, the pair potential reads: 
√
(2) (2) (2) (2)
VQ(Ωˆi, Ωˆj)= −E Δ(Ωˆij)+ 2λ[Δ(Ωˆij)+Δ(Ωˆij)] + 2λ2Δ(Ωˆij) , (3.4)
00 2002 22 
which is the parametrization of the model as presented in [64] and [104]. Clearly, the minimum of (3.1) is achieved when all of the respective molecular axes of Q(Ωˆi) and Q(Ωˆj) coincide, i.e. when: 
ˆ
aˆi aˆj,bi ˆbj,cˆi cˆj. (3.5) 
Thus, the ground state of (3.1) is the biaxial nematic phase NB (Fig. 3.1b). 

ˆˆ
a) exemplary orientations of the molecular tripods Ωˆi = {aˆi,bi,cˆi} and Ωˆj = {aˆj,bj,cˆj}with respect to the laboratory frame of reference {ˆy, ˆ
x, ˆz}; b) a ground state of the potential (3.1), in which aˆi aˆj, ˆbi ˆbj,cˆi cˆj. The molecular axes are scaled by the magnitudes of the corresponding eigenvalues of the molecular quadrupolar tensors (a red arrow denotes a negative eigenvalue). The pictured quadrupolar tensors are for λ =0.3. 
3.2.1.1 Properties of Q(ˆ
Ω) ˆ
Consider a molecular tensor Q(Ω) ˆconstructed in the molecular basis Ωˆ= {ˆb, ˆ
a, c} for an arbitrary particle. Since this tensor is traceless, we cannot directly equate its eigenvalues with spatial dimensions of a tangible object. To establish an intuitive framework, from now on we adopt a convention that the sign of the largest-modulus non-degenerate eigenvalue of Q(ˆ
Ω) specifies the molecular long axis (+) or the largest face (−). The axis corresponding to this eigenvalue is considered the primary molecular axis, while the second largest-modulus eigenvalue (if non-degenerate) identifies the secondary molecular axis. In the diagonal form, Q(ˆ
Ω) reads: 
⎞
⎛ 
−√1+ λ 00 
6 
0 −√1
⎜⎜⎝ 

⎟⎟⎠

R

{ˆ
a,
ˆc},
b,ˆ
(3.6)

Q(ˆ
Ω) = RT 
{ˆ
a,
ˆ
− λ 

0 

b,cˆ} 
6 
2
00 
3 
where Rˆis a rotation matrix, carrying the laboratory frame to the molecular frame 
{a,ˆb,cˆ}
of reference. From the above form we can now see that in the limit cases of the pa
3
rameter λ ∈ [0, ] a uniaxial tensor is recovered. Namely, for λ =0 the eigenvalues 
2 

Ω) and its eigenvalues in the di
agonal basis for a) λ<

1√ 6
, b) λ>

1√ 6 
and c) λ 

=

1√ 6 
. 

The vectors are scaled by 

magnitudes of corresponding eigenvalues. A red arrow means that the corresponding eigenvalue is negative. Surface described by r = exp (Q · uˆ⊗ uˆ) where uˆis a unit vector in polar coordinates. 
2 
3
ˆ
}, so that Q(ˆa,
Ω) is degenerate in the {ˆb}-plane and cylindrically 
c. The largest-modulus non-degenerate eigenvalue corresponds to cˆand 
1√
1√
are {−symmetric along ˆ
, −

,

6
6
is positive. Therefore, this case corresponds to the prolate molecule with primary axis cˆ. (See Fig. 3.2a.) 
Conversely, if λ = 

3 2
, the eigenvalues are {

2 3
, −2

22
}, axes {ˆ

a, cˆ} are degenerate 

,

33
and the largest-modulus non-degenerate eigenvalue is negative and corresponds to ˆb. Thus, the corresponding molecular shape is a flat disk with primary axis ˆb. (See Fig. 3.2b.) 
In the self-dual point λ 

= 

degenerate, the axes ˆb and cˆare equivalent. The corresponding object has both a pro
1√ 6
, the eigenvalues are {0, −

22
} and, while non
,

33
found long axis, identified by cˆ, and a large face, identified by ˆb. This case corresponds 
√ 
to biaxial bricks of dimensions (L, LW,W ), as identified by Straley. (See Fig. 3.2c.) 
3.2.1.2 Duality transformation 
Thedualitytransformationfor model(3.1) isperformed by simultaneouslyrescaling both λ and thedimensionless temperature t ≡ kBT/E, while cyclically permuting the axessuch 
ˆˆ
b, ˆc, ˆb} [108]:
that {ˆ
a, 
c}→{ˆ

a, 

λ → 

3 
2
− λ √ 
1+ 6λ (3.7) 4t 
t →√ . 
(1+ 6λ)2 
Such transformation leaves the Boltzmann factor e−H/t invariant. (3.7) is identity when 
λ = 

1√ 6 
. 

3.2.1.3 Phases and phase transitions 
The isotropic-uniaxial nematic phase transition is of first order, occurring through spontaneous breaking of the global rotational symmetry of the Hamiltonian (3.1), as in the Maier-Saupe and Lebwohl-Lasher models. In the uniaxial nematic phase the long-range ordering is of the molecular axis corresponding to the largest-modulus eigenvalue of 
Q(ˆ
Ω). Therefore, a nematic prolate NU+ phase is observed for λ<
1√ 6 
and a nematic 

oblate NU− is observed for λ>
1√ 6
, with the primary uniaxial director nˆdue to the 

long-range ordering of the cˆand ˆb axes respectively. For these cases, as the temperature 
is lowered in the nematic phase, a second-order phase transition to the biaxial nematic 

NB 
phase occurs, through spontaneous breaking of uniaxial symmetry. At λ = 

1√ 6 
nei
ther of cˆand ˆb is privileged and ordering of both is simultaneous, hence there is a direct transition from the isotropic (I) to the biaxial nematic phase (NB), known to be second order [104, 105]. The phase diagram for the model (3.1) is presented in Fig. 3.3. 

1
correspond to mean-field predictions. At the self-dual point λ = √ ≈ 0.408248 there 
6 
exists a direct isotropic-biaxial nematic second order transition. Here the symbols I, N+, N− and B stand for the isotropic, nematic uniaxial prolate, nematic uniaxial oblate and nematic biaxial phases, throughout the present text denoted as I, NU+, NU− and NB respectively. Reprinted FIG. 2 with permission from F. Biscarini, C. Chiccoli, P. Pasini, F. Semeria, and C. Zannoni. Physical Review Letters, 75 (9):1803–1806, 1995. 
Copyright 1995 by the American Physical Society.) 


3.3 Dispersion model with external field 
In nematic liquid crystals application of an external electric or magnetic field affects the orientation of the sample [2, cc. 3.2-3.3]. Most nematogenic molecules are diamagnetic, as they in general contain one or more benzene rings (e.g. MBBA, see Fig. 3.4), which respond to the application of an external magnetic field by an induced repellent eddy current, if the field is at an angle to the benzene ring plane. Therefore, in the case of magnetic fields, the effect works generally to adjust the director n 
to align in parallel to the external magnetic field Hn. 

The magnetization induced in a nematic sample by an external magnetic field Hnreads: 
nn
M = χ⊥H +(χ − χ⊥)(Hn· n
)nn. (3.8) 
χ and χ⊥ are the anisotropic components of the bulk diamagnetic polarizability tensor, corresponding to the long molecular axis and the axis perpendicular to benzene ring plane. The difference Δχ = χ − χ⊥ is the diamagnetic anisotropy of the sample. The corresponding contribution to energy is given by: 
� H 
nn
ΔEM = − M · dHn= − 1 χ⊥H2 − 1Δχ(n 
· H)2 . (3.9)
22
0 
The first term is independent of n
, so it can be discarded. The second term promotes nn
n
H if Δχ> 0 and n 
⊥ H if Δχ< 0. For uniaxial nematics the alignment tensor is Q = S(n 
⊗ n 
− 1 I), so the second term can be rewritten:
3 
    
− 1 2 Δχ(n 
· nH)2  = − 1 2 Δχ HiHjQij +  1 3 H2 ,  (3.10)  

where i, j = {x, y, z}. Again, the term proportional to H2 can be dropped, and we have: 
1 
ΔE ' = − ΔχHiHjQij. (3.11)
M 
2 
A similar term to (3.11) appears when the electric field is considered instead. 

3.3.1 Model 
A similar discussion to (3.8-3.11) can be conducted at the microscopic level. The contri<BR>bution to the potential energy of particle i then reads [70]: 
ΔV (Ωˆi)= −EΔEhαhβQ(ˆ(3.12)
Ωi)αβ. 
We distinguish the general vector field nh because of the universal character of (3.12), which can also express interaction with an electric field (to an extent where piezoelectric effects can be neglected [110]). The analogue of bulk diamagnetic anisotropy is the anisotropy of the molecular (dielectric or diamagnetic) field polarizability, EΔE. In our discussion we distinguish only the sign of the anisotropy of the field polarizability and let ΔE = ±1, either assuming that the anisotropy is equal to the quadrupolar coupling constant E, or that it renormalizes h2. Note that we use the full biaxial molecular tensor 
Q(ˆ
Ω), as defined in (3.2), and assume a proportionality for the molecular polarizability tensor for particle i: χi (3.13)
αβ ∝ EΔEQ(Ωˆi)αβ, 
which is a model assumption. Together with (3.1), the dispersion model for biaxial nematics in the presence of an external vector field is given by the Hamiltonian: 
⎛  ⎞  
H = −E ⎝1 2  N  Q(ˆΩi) · Q(ˆΩj) + ΔE  N  hαhβQ(ˆΩi)αβ⎠ .  (3.14)  
(ij)  i  

We take the first sum over the nearest neighbors, i.e. for every i we take only z =6 neighboring particles on a simple cubic lattice with periodic boundary conditions. This model corresponds to the general model (2.11) for the case when the parameters λ, E and h are non-zero. In our discussion we will assume that the model (3.14) is studied with respect to a dimensionless temperature t ≡ kBT/E (which is equivalent to putting E =1). 

3.3.2 Field-induced effects 
While the model (3.1) has global rotational symmetry, in the model with non-zero exter<BR>nal field (3.14) the rotational symmetry is broken by definition. Therefore, the Hamil<BR>tonian (3.14) has global symmetry with respect to rotation around the vector nh, which determines the preferred alignment of one of the molecular axes and, consequently, one of the directors. Considering a sample particle for which the orientation is given by 
ˆˆ
Ω= a, c}, there are two main microscopic effects, consistent with the minima of 
{ˆb, ˆ(3.12), when h �=0, namely: 
• 	
for ΔE = +1 the molecular cˆaxis (which corresponds to the largest, with respect to the modulus, positive eigenvalue of Q(ˆh

Ω), regardless of λ) aligns parallel to n(Figs. 3.5a and c); 

• 	
for ΔE = −1 the molecular ˆb axis (which corresponds to the largest, with respect to the modulus, negative eigenvalue of Q(ˆh


Ω), regardless of λ) aligns parallel to n(Figs. 3.5b and d). 
What precisely happens to the orientation of the director basis depends on λ and the considered zero-field phase, to which a non-zero field is introduced. We identify two classes of macroscopic effects, distinguished by how the field affects the director of the zero-field uniaxial nematic phase. Namely, considering the above discussion of the orientation of molecular axes, the field can make the primary uniaxial director nˆeither parallel or perpendicular to nh. In the first case, the field reinforces uniaxiality, while in the second case biaxiality is induced, because the field forces rotation of the uniaxial director nˆto align perpendicularly to nh and orients the secondary molecular axis (and 
n
therefore the secondary director mˆ) in parallel to h. 
The uniaxiality-promoting effect occurs for λ< 

√
√
1166 
∧ ΔE = +1 or λ> 

∧ ΔE = −1 

(Figs. 3.5a and b). <BR><BR>
The biaxiality-promoting effect occurs for λ< 

√
√
1166 
∧ ΔE = −1 or λ> 

∧ 

ΔE = +1 

(Figs. 3.5c and d). 
3.3.2.1 Uniaxiality-promoting effects 
Cases of λ< 

√
16 
∧ ΔE =+1 (Fig. 3.5a) and λ> 

√
16 
∧ ΔE = −1. (Fig. 3.5b). 

A non-zero field aligns the director nˆparallel to the field in all phases for this particular case. Firstly, in the isotropic phase the field induces uniaxial ordering. The increase in uniaxiality (e.g. in the uniaxial order parameter S, as defined in (1.2)) is proportional to h2 and the resulting phase is called paranematic. Depending on λ this uniaxial phase 
can be either prolate (NP +) for λ< 
1√ 
6 
, or oblate (NP −) for λ> 
1√ 
6 
. 

When considering thezero-fielduniaxial nematicphase, the uniaxial ordering is increased by the application of the field. The NP − NU phase transition is of first order. It is known that the transition temperature rises with increased field and it is predicted by phenomenologicaltheoriesthatthetransitionshouldterminate atacriticalpoint[70,96]. 
In the biaxial nematic phase the field does not reduce symmetry, although Landau-de Gennes theory predicts that the NU − NB transition temperature may either increase or decrease [103]. 

a) λ<

b) λ>

c) λ<

d) λ>

∧ ΔE =+1. 6
6
6 
6 
1√ 1√ 1√ 1√
∧ ΔE = +1, 

∧ ΔE = −1, 

∧ ΔE = −1, 

Thearrowsrepresenteigenvectorsof Q(ˆ
Ω) andarescaledbythemagnitudesofrespective eigenvalues. A red arrow means that the corresponding eigenvalue is negative. The dashed arc illustrates the degeneracy of the minimum in the plane perpendicular to nh. Surface described by r = exp (Q · uˆ⊗ uˆ), where uˆis a unit vector in polar coordinates. 
3.3.2.2 Biaxiality-promoting effects 
Cases of λ< 

√
16 
∧ ΔE = −1 (Fig. 3.5c) and λ> 

√
16 
∧ ΔE =+1. (Fig. 3.5d). 

In the isotropic phase the field orients the director nˆparallel to the field, as in the 

previous case. Depending on λ, the paranematic prolate phase (NP +) occurs for λ>
1√ 6 
and the paranematic oblate phase (NP −) is found for λ<
1√ 6 
. 

n
In the uniaxial nematic phase the field forces the director nˆto align in perpendicular toh and biaxiality is induced, with the secondary director mˆparallel tonh. Upon the NP −NB 
n
transition, the director nˆis reoriented from being parallel to h to being perpendicular to hˆ. This transition is initially first-order, but is predicted by phenomenological theories to experience a tricritical point for high enough fields and then to continue as second-order [70, 96], when the field is sufficiently large so that field-ordering becomes dominant. 
The biaxial nematic phase, which exists for no field, blends seamlessly with the biaxial nematicphasewhichisinducedinplaceoftheuniaxialnematic. The NU −NB transition, being second order, simply vanishes as the field is introduced, leaving a smooth increase in biaxial ordering as temperature is lowered. 
3.3.2.3 Field effects for λ = √1
6 
Theself-dualpointisaspecialcasebecauseofthehighsymmetryinvolved. For h =0 the molecular axes cˆand ˆb are equivalent and there is no director along which the uniaxial ordering is predominant, as it can only be field-induced. For h �c is preferred to align 
=0ˆ
nn
in parallel toh for ΔE = +1, while ˆb is preferred to align in parallel toh for ΔE = −1. In the isotropic phase, this induces a paranematic prolate (NP +) phase and a paranematic oblate (NP −) phase in the respective cases. In the self-dual case the biaxial nematic phase for h =0 is a phase with maximum biaxiality (that is, the alignment tensor has no uniaxial part), in which the induced uniaxial ordering for h �
=0 can actually reduce the level of biaxial ordering, thus lowering the NP − NB transition temperature. 

3.3.3 Duality transformation 
Thediscussedfield-inducedeffectshintuponasymmetryofthemodel(3.14)withrespect to the permutation of molecular axes and a change of λ from the prolate to the oblate case, coupled with a switch in the sign of ΔE. Such self-dual transformation for (3.14) is performed by extending on the transformation (3.7) for model (3.1) by simultaneously rescaling also the magnitude of the field. Namely, the argument of the Boltzmann factor 
ˆˆ
e−H/t is invariant when the axes are cyclically permuted such that {ˆb, ˆc, ˆb}
a, c}→{ˆa, and dimensionless temperature t = kBT/E, λ and ΔEh2 are rescaled: 
λ →  3 2 − λ √  
1 +  6λ  
t →  4t (1 + √ 6λ)2 .  (3.15)  
ΔEh2 →  −2ΔEh2 √  
1 +  6λ  



3.4 Tools and methods 
3.4.1 Order parameters 
3.4.1.1 Alignment tensor 
The description of order in biaxial nematics is primarily achieved by investigating the properties of the alignment tensor, which, similarly to the molecular quadrupolar moment, is composed of uniaxial and biaxial contributions: 
√
(2) (2)
Q ≡ T(Ωˆ' )+ 2λT(Ωˆ' ), (3.16)
02 
where the overline denotes the ensemble average and Ωˆ' is, as of yet, an unspecified cyclic ˆ
permutation of the molecular basis {a,ˆb, cˆ}, for which there are three possibilities. The 
(2) 
Ω ' (2.8-2.9). The proper choice of Ωˆ', for which a consistent Q is obtained, is made such that tensors Tm (Ωˆ' ),m =0, 2 are constructed by assuming {eˆx,eˆy,eˆz}≡ ˆin the definitions 
(2) (2)
theaverageuniaxialtensor T(Ωˆ' ) andtheaveragebiaxialtensor T(Ωˆ' ) unequivocally 
02 
determine the primary uniaxial director nˆand the two secondary (biaxial) directors mˆand ˆl respectively, consistently with the qualitative predictions made in the previous n
sections, depending on the field h and parameters λ and ΔE: 
1
• For λ �√ theappropriateuniaxialtensoristheoneconstructedupontheprimary 
= 
6 
molecular axis (while the biaxial tensor is built out of the remaining two), thus 
11
Ωˆ' = {ˆb, ˆ√ and ˆ= {ˆa, ˆ.
a, ˆc} for λ< Ω ' c, ˆb} for λ> √ 
66 
1
• For λ = √ and h =0 in the biaxial phase the only viable choice is Ωˆ' = {ˆb, c,ˆaˆ},
6 
for which the biaxial tensor is constructed using the (equivalent) axes cˆand ˆb. In the other two options one of them is being distinguished completely artificially. 
(2)
The average uniaxial tensor T(Ωˆ' = {ˆb, c,ˆaˆ}) is purely uniaxial and is non-zero 
0 
only because the ordering of cˆand ˆb necessarily causes ordering of aˆ. 
1
• For λ = √ and h =0 �the field prefers rotation of one of the molecular axes parallel 
6 
to nh, i.e. cˆfor ΔE = +1 and ˆb for ΔE = −1. In this case the average of T(2)(Ωˆ' )
0 
Ω '
with the choice of ˆ= {ˆb, c,ˆaˆ}, as for h =0, is not consistent with the uniaxial director and also does not account for the increase in uniaxial ordering along nh. 
Ω ' ˆΩ ' ˆ
Therefore, ˆ= {ˆb, ˆ> 0 and ˆc, ˆb} for ΔEh2
a, c} for ΔEh2 = {ˆa, < 0. 
In conclusion, the alignment tensor which is consistent with the qualitative predictions is defined as: 
Q ≡ 

⎧ ⎪
⎪
⎨ ⎪
⎪
⎪⎩ 

√
(2) (2) 	11
ˆ	ˆ
T({ˆb, ˆ2λTa, c}) for λ< √ and λ = √∧ ΔEh2 
a, c})+ ({ˆb, ˆ	> 0
0	2 
66 
√
(2) (2) 	11
ˆ	ˆ
T({ˆa, b})+ 2λTc, ˆb}) for λ> √ and λ = √∧ ΔEh2 
c, ˆ({ˆa, 	< 0
0	2 
66 
. 

√ 
1
(2) 	(2)
T({ˆb, c,ˆaˆ})+ 2λT({ˆb, c,ˆaˆ}) for λ = √∧ h =0
0	2 
6 
(3.17) 
By following the definition (3.17), we can calculate the director frame of reference {ˆl, ˆn} as the eigenbasis of Q, identifying the directors by the moduli of its respec
m, ˆ
tive eigenvalues: {ˆl, ˆn} : {|ωˆlm|≤|ωnˆ|}. (3.18)m, ˆ|≤|ω ˆ
Because in (3.17) the averaged tensors are traceless, also Tr(Q)= n + ω ˆ=0.
ωˆm + ωˆl The properties of the phase in which Q is measured can be inferred from the eigenvalues: 
• In the isotropic (I) phase ωnˆm =0 (no directors can be identified).
= ω ˆ= ωˆl 
• 	In the uniaxial nematic or paranematic prolate phase (NU+ or NP +), ωnˆ> 0,ω ˆ= 
m 
ωˆ= −ωˆ/2 (nˆis well-defined, while Q is degenerate in the {ˆl, mˆ} plane).
ln
• In the uniaxial nematic or paranematic oblate phase (NU− or NP −), ωnˆm
< 0,ω ˆ= ωˆ= −ωnˆ/2 (nˆis well-defined, while Q is degenerate in the {ˆl, mˆ} plane).
l 
1
• 	In the biaxial nematic phase for (NB) for λ �√ , ωˆ�ω ˆ,ωˆ= −(ωˆ+ ω ˆ)
= n = mnm
l
6 1
(all three directors are well-defined). For the isolated case of λ = √ ∧ h =0,
6 ωnˆml 
= −ω ˆ,ωˆ=0 (directors nˆand mˆare interchangeable). 
• 	In addition, we can use the invariant parameter w, as defined in (1.21), which takes different values in the various phases: w = +1 in prolate phases (NU+ or NP +), w = −1 in oblate phases (NU− or NP −) and |w| < 1 in the biaxial nematic phase (NB). In particular w =0 in the maximally biaxial phase. In the isotropic phase w is undefined. 
3.4.1.2 Definition of order parameters 
The natural order parameters are the coefficients for decomposing Q in the basis of symmetrized irreducible tensors in the director basis: 
31
(2)
T({ˆl, ˆn})= n ⊗ nˆ− I)
0 m, ˆ(ˆ23 (3.19)
1
(2)
T({ˆl, ˆn})= √ (ˆl ⊗ ˆm ⊗ mˆ).
2 m, ˆl − ˆ
2
The decomposition is constructed using the D2h-symmetrized identity operator (2.7) for l =2 (for which case the seniority index is irrelevant), constructed in the director basis: 
Q =  m=0,2 Q · T(2) m ({ˆl, ˆm, ˆn}) ⊗ T(2) m ({ˆl, ˆm, ˆn}),  (3.20)  
where the coefficients:  
q0 q2  = =  Q · T(2) 0 ({ˆl, ˆm, ˆn}) Q · T(2) 2 ({ˆl, ˆm, ˆn})  (3.21)  

are the uniaxial (q0) and biaxial (q2) order parameters, normalized such that |q0|∈ [0, 1]√ 
and |q2|∈ [0, 2λ]. With these parameters we can express the eigenvalues of Q, i.e.: 
2 q0q2q0q2
{ωˆ= q0,ω ˆ= −√ −√ ,ωˆ= −√ + √}. (3.22)
nm 	l
362 62
From the eigenvalues of Q (3.22) and the discussion from the previous section we infer the following properties of the order parameters: 
• 	
In the isotropic (I) phase q0 = q2 =0. 

• 	
Intheuniaxialnematicorparanematicprolatephase(NU+ orNP +), q0 > 0,q2 =0. 

• 	
In the uniaxial nematic or paranematic oblate phase (NU− or NP −), q0 < 0,q2 =0. 

• 	
In the biaxial nematic phase (NB) either q0 > 0,q2 < 0 or q0 < 0,q2 > 0, which is enforced by Tr(Q)=0. 


3.4.1.3 Invariant parameter w 
In the cases where the increases in ordering are very small, as in the field-induced NP phases, the invariant parameter w (1.21), as introduced in section 1.6, provides unequiv<BR>ocal identification of the phase, assuming a value of ±1 for arbitrarily small uniaxial 
ordering. Expressed through order parameters (3.21), the parameter reads: <BR><BR>
32
q0 − 3q0q
2 
w = . 	(3.23)
22)3/2
(q0 + q
2

3.4.2 	Mean-field method 
Byapplyingtheformalismrecollectedinsection(1.7)tothepotentialofthemodel(3.14), we write the mean-field equations, which enable us to calculate the order parameters q0 and q2. The effective potential is obtained as follows: 
z=6 Veff (ˆΩi)Q(ˆΩ) Ω),
Ω) = −E Tr Peq(ˆΩi) · Q(ˆ= −EzQ · Q(ˆ(3.24)
ˆ
Ωi
i 
ˆ
where z is the coordination number (number of nearest neighbors) and Ωˆ{ˆb, cˆ} is
= a, 
the orientation of the molecular basis for a sample particle, expressed with respect to the laboratory frame of reference. By applying the decomposition of the alignment tensor 
(3.20) to (3.24), we obtain: 
(2) 	(2)
Veff (ˆ02 l, ˆ	(3.25)
Ω) = −Ez[q0T({ˆl, m,ˆnˆ})+ q2T({ˆm, nˆ})] · Q(Ω)ˆ. 
The above definition can be explicitly expressed in terms of symmetrized Wigner functions, discussed in section 2.1: 
(2) 	(2)
ΔΩ)= T(2) m, ˆm ' (ˆ	(3.26)
(ˆ({ˆl, ˆn}) · TΩ),
mm ' m 
where we use the definition (2.6), substituting the director basis for Ωˆ1 ≡{ˆm, ˆ
l, ˆn} and ˆ
the molecular basis for a sample particle for Ωˆ2 ≡{ˆb, ˆ
a, c}. Using (3.26) in (3.25), while observing the form for the tensor Q(ˆ
Ω) (3.2), we obtain: 
√	√
(2) (2) (2) (2)
Ω) = −Ez ΔΩ)+ 2λΔ(ˆΔΩ)+ 2λΔ(ˆ.
Veff (ˆq0 00 (ˆ02 Ω) + q2 20 (ˆ22 Ω) (3.27) 
The complete mean-field Hamiltonian and equilibrium free energy per particle read: 
HMF (ˆ= Ω) − EΔEhαhβQ(ˆ	(3.28)
Ω) 	Vef (ˆΩ)αβ, z 
22
f = F/N =(q0 + q2) − t log Z. 	(3.29)
2
Themean-field partition function and equilibrium particledistribution function are given by: 
Ωe−βHMF (ˆ
Z = dˆΩ), (3.30) 
Ω) = e Ω), (3.31)
Peq(ˆZ−1 −βHMF (ˆ
where β =(kBT )−1. The order parameters can be calculated as mean-field averages of combinations of symmetrized Wigner functions, as given below: 
√
(2) (2)
dˆ(ˆ(ˆ(ˆ
q0 = Ω) ΔΩ)+ 2λΔΩ) ,
ΩPeq00 02 
� √ (3.32) 
(2) (2)
q2 = dΩˆPeqΩ) (ˆ2λΔΩ) .
(ˆΔΩ)+ (ˆ
20 22 
q0 and q2 are obtained using two methods, described in the paragraphs that follow. 
3.4.2.1 Self-consistent approach 
In the first procedure the equations (3.30-3.32) are solved iteratively. We parametrize Ωˆwith Euler angles and the integration is carried out using a Monte Carlo method with importance sampling, through the Vegas algorithm from the Cuba Library, version 2.1 [111]. The initial values of parameters q0 and q2 for a given temperature t are supplied randomly with a uniform distribution on [−1, 1], and subsequently we observe which seed yields convergence to the lowest free energy per particle f (3.29). The equations are iterated until the variance of the last 10 iterations drops below 10−10. It typically takes 20 ∼ 30 iterations for convergence when t is far from the transition temperature, while as much as 100 ∼ 200 iterations are needed in the phase transition vicinity. Because the integration is approximate, the behavior of order parameters is correctly recovered only away from transition temperature, while in the vicinity of the phase transition the order parameters are smoothed out, which makes this method impractical for distinguishing between transitions of different order. However, the self-consistent method is very useful in studying the NU −NB transition, for which the order is known. This transition occurs at relatively low temperatures, where the high-temperature expansion method, discussed below, proves difficult to apply. With the self-consistent method we successfully perform calculations for (dimensionless) temperatures as low as 0.1. It should be noted, however, that anything that is obtained in such low temperatures is only a property of the model, 
|t−tINU |
since real nematic phases are stable for higher temperatures, for which tINU � 0.2, where tINU is the temperature of the isotropic to nematic uniaxial phase transition. 
3.4.2.2 High-temperature expansion 
We expand the mean-field partition function (3.30) around βE = t−1 =0 (t being the dimensionless temperature) in the expression for the mean-field free energy per particle 
(3.29):<BR><BR>
 
 
n 
i
−i 1
z 
22 dˆ
(q0 + q2) − t log Ω
Ω) . (3.33)
−HMF (ˆ
fn 
= 

t

2

i! 
i=0 
The procedure is to choose a suitable expansion order n and integrate term-by-term, leaving q0 and q2 as free parameters. The terms are integrable polynomials in trigonometric functions of Euler angles. q0 and q2 are found for a given t by minimization of the resulting approximate free energy per particle fn: 
⎧ ⎪
⎪ 
⎪
⎪
⎨ ⎪ 
⎪
⎪
⎪⎩ 

∂ ∂q0 fn =  0  
∂ ∂q2 fn = ∂2 ∂q2 0 fn >  0 0  (3.34)  
∂2 ∂q2 2 fn >  0.  

We use n =6 through n = 14 as needed. Because of the factor t−i, the series in (3.33) converges very slowly in low temperatures, where higher terms become important and even for orders as high as n = 14 the approximation breaks down at t ≈ 0.8, which (save 
1
for the vicinity of λ = √ ) is usually much higher than the NU − NB transition. On the 
6 
other hand, in t> 1, as few as n =6 terms provide a very good approximation, therefore the expansion (3.33) proves very useful in studying the I − NU transition. Because the leading terms are exact and not approximated, as in the self-consistent method from the previoussection, theexpansionmethodcorrectlyrecoversorderparameterdiscontinuities at the first-order phase transition. 

3.4.3 Monte Carlo simulations 
We assume a simple cubic lattice of N = 16 × 16 × 16 particles with periodic boundary conditions. Each molecule interacts only with the z =6 nearest neighbors. For each separate temperature the simulation runs independently, starting from a completely biaxial initial configuration of the lattice, such that the molecules are aligned with their molecular axes parallel to the laboratory frame of reference {ˆy, ˆ
x, ˆz}, which is parallel to the lattice edges. That is, in the initial state, for each lattice site i: 
aˆi ˆˆbi ˆci ˆ(3.35)
x, y, ˆz. 
Furthermore, we put the field parallel to the zˆaxis: 
n
h ˆ(3.36)
z. 
Subsequently, the Monte Carlo step radius r, defined in (1.49), is self-consistently ad<BR>justed so that the average number of accepted trial moves is between 40% and 50%. As the system evolves, the acceptance may increase or lower. In particular, in the isotropic phase it can rise to as much as 60% when the radius is not adjusted. We decide to keep the radius constant during the simulation, as adjustment during simulation could lead to violation of detailed balance [112]. 
One Monte Carlo lattice sweep (which, generally, is called a Monte Carlo cycle in literature) proceeds by performing trial moves on N randomly chosen molecules. For a molecule i the change in energy upon a trial rotation Ωˆ' is calculated as: 
i 
z=6
1 
ΔEi = V (Ωˆi, Ωˆj) − V (Ωˆ' i, Ωˆj)+ΔV (Ωˆi) − ΔV (Ωˆ' i) , (3.37)
2 
j(i) 
where the two particle potential and external field interaction are as defined in (3.1) and (3.12). The trial move is accepted according to the Metropolis algorithm discussed in section 1.8.1. 
As many as 100k cycles are needed for the thermalization process. For calculating averages, or production, we perform from 100k to 240k cycles. As previously discussed in section 1.8, one needs to discard correlated states when calculating averages. We find that taking every tenth equilibrium configuration is a safe choice, so in the end we consider as many as 10k to 24k configurations when calculating thermodynamic quantities. Thethermodynamic average ofageneral microscopicproperty Ai, definedfor each lattice site i =1,...,N, is calculated as: 
A = 1 M  M  1 N  N  Ai  = 1 M  M  ˜Ak,  (3.38)  
k=1  i=1  k  k=1  

where M is the size of the Monte Carlo statistical sample, in which k numbers the elements and A˜k is the k-th instantaneous lattice average. Furthermore, we define the 
Monte Carlo fluctuation of Ai as: 
2
MM
11 
χA = A˜2 M − A˜2 = A˜k 2 − A˜k , (3.39)
M 
MM 
kk 
where it is presented that the average · M is calculated as A˜M = 1  M A˜k.
Mk 
From Monte Carlo simulations we obtain the alignment tensor Q, as defined in (3.17), from which the director basis {ˆl, ˆn} is calculated by diagonalization and labeling of 
m, ˆthe respective eigenvectors according to (3.18). Order parameters q0 and q2 are then calculated from the definition (3.21). 
3.4.3.1 Specific heat and susceptibilities 
To aid indication of phase transitions, we calculate the specific heat, which has an important property that it is singular at a phase transition. In the canonical ensemble the specific heat at constant volume is proportional to the fluctuation of energy: 
�  �  
CV  =  ∂ E ∂T  V  = kB T 2 ( E2  − E 2).  (3.40)  

We express CV through fluctuation of the per-particle energy: 
z=6
1 
ˆ
εi = V (Ωˆi, Ωj)+ΔV (Ωˆi). (3.41)
2 
j(i) 
The fluctuation of per-particle energy is calculated as in definition (3.39): 
  
c = Nε˜2 M − ε˜2 , (3.42)
M
where ε˜is the lattice average of εi. CV can be calculated from c: 
kB
CV = N c. (3.43)
T 2 
For better clarity of presentation, we will use c rather than CV , because the factor 1/T 2 makes the peaks for the low-temperature transitions less visible. 
It is also useful to consider susceptibilities of order parameters, which too are singular at a phase transition. Statistical physics relates the susceptibilities to fluctuations of order parameters. Since we do not compute the instantaneous director basis for each equilibrium state, instead of calculating fluctuations of q0 and q2, we calculate the fluctuations of scalar order parameters, given by instantaneous lattice averages of molecular uniaxial 
˜(2) ˜(2)
and biaxial tensors T0 and T2 : 
˜(2) ˜(2) ˜(2) ˜(2) 2
χ0 = T· T− T· T, (3.44)
00 00 
˜(2) ˜(2) ˜(2) ˜(2) 2
χ2 = T· T− T· T. (3.45)
22 22 
˜(2) ˜(2)
Quantities(3.44-3.45)areproportionaltothesusceptibilitiesoforderparameters Tm · Tm , 
but we will refer to them simply as “susceptibilities of uniaxial and biaxial order parameters”, as we do not make use of their quantitative interpretation. 


3.5 Results for the Luckhurst-Romano model 
As a test of our methods, we present the results recovered for the original Luckhurst-Romano model (3.1), i.e. when h =0 in (3.14) for three representative values of the parameter λ. The results are in agreement with previous research [104, 105]. 
The case of λ =0.3 corresponds to prolate biaxial molecules. Results are presented in Figs. 3.6a-c. The transition temperatures are given in Table 3.1. The phase transition sequence with lowering temperature is as follows: 
• 	
In the high-temperature regime the isotropic phase (I) is present. In this phase the order parameters q0 and q2 are zero and the invariant parameter w is undefined (because of Tr(Q2)=0). 

• 	
At t = tINU+ a weakly first-order phase transition occurs from the isotropic (I) to the uniaxial nematic prolate (NU+) phase. The energy fluctuation c and uniaxial order parameter susceptibility χ0 experience sharp peaks at the transition temperature, while the uniaxial order parameter q0 begins to rise and w approaches 1. The biaxial order parameter q2 is zero. tINU+ ≈ 1.08 in Monte Carlo, tINU+ ≈ 1.48 in mean-field. 

• 	
At t = tNU+NB a second-order phase transition occurs from the uniaxial nematic prolate (NU+) to the biaxial nematic (NB) phase. The energy fluctuation c and the biaxial order parameter susceptibility χ2 experience peaks at the transition temperature, while the biaxial order parameter q2 begins to rise (modulus-wise), whereas w begins to decrease from 1. tNU+NB ≈ 0.43 in Monte Carlo, tNU+NB ≈ 


0.59 in mean-field. 
The case of λ =0.5 corresponds to oblate biaxial molecules. Results are presented in Fig. 3.7a-c. The phase transition sequence with lowering temperature is as follows: 
• 	
As in the previous case, in the high-temperature regime the isotropic phase (I) is present. In this phase the order parameters q0 and q2 are zero and the invariant parameter w is undefined (because of Tr(Q2)=0). 

• 	
At t = tINU− a weakly first-order phase transition occurs from the isotropic (I) to the uniaxial nematic oblate (NU−) phase. The energy fluctuation c and uniaxial order parameter susceptibility χ0 experience sharp peaks at the transition temperature, while the uniaxial order parameter q0 begins to rise (modulus-wise) and w approaches −1. The biaxial order parameter q2 is zero. tINU− ≈ 1.35 in Monte Carlo, tINU− ≈ 1.85 in mean-field. 


q 0 , q 2 
c 

0.2 

2 
0 
-0.5 
1
-0.2 
-0.4 
-1 
0 
t 
c ab 
5
1 
1 
0.8 
4 
0.5 
0.6 
3
0.4 

w

0 

Figure 3.6: Results for λ =0.3 and h =0. Monte Carlo transition temperatures are marked with vertical dashed lines, while the phases are labeled above the plot. a) Monte Carlo (points) and mean-field (lines) uniaxial and biaxial order parameters q0 and q2 (3.21). Note how the change in mean-field q0 is abrupt, indicating the first-order character of the I − NU+ transition. b) Monte Carlo (points) and mean-field (lines) invariant parameter w (1.21) and Monte Carlo energy fluctuation c (3.42) (in green). When calculating w from Monte Carlo simulations results for which Tr(Q2) < 10−3 are discarded, which happens in the isotropic phase. c) Monte Carlo fluctuations of uniaxial (χ0) and biaxial (χ2) order parameters. In Monte Carlo results, at temperatures directly above tINU+ ≈ 1.08, minute increase in uniaxiality and biaxiality is visible (see w, q0 and q2), which is identified as a finite-size effect, often observed for such system sizes [104, 105]. The mean-field lines are composed of the high-temperature expansion result for n =8 (see section 3.4.2.2) for t ≥ 1 and of the self-consistent result (see section 3.4.2.1) for t< 1. 
• 	At t = tNU−NB a second-order phase transition occurs from the uniaxial nematic oblate (NU−) to the biaxial nematic (NB) phase. The energy fluctuation c and the biaxial order parameter susceptibility χ2 experience peaks at the transition temperature, while the biaxial order parameter q2 begins to rise, whereas w begins to increase from −1. tNU−NB ≈ 0.62 in Monte Carlo, tNU−NB ≈ 0.86 in mean-field. 
1
In the self-dual case λ = √ there exists one phase transition, from the isotropic (I) to 
6 
the biaxial nematic (NB) phase, which is of second order. Monte Carlo and mean-field results are presented in Fig. 3.8. The qualitative difference of the I −NB transition from the first-order I − NU and second-order NU − NB phase transitions comes from the fact 
1
that for λ = √ two directions of symmetry breaking (two perpendicular directors) are 
6 
established simultaneously. In Monte Carlo results this leads to an apparent ambiguity of the transition temperature. In fact, the fluctuation-measuring quantities c, χ0 and χ2 are peaked at a temperature which is ∼ 18% lower than the temperature at which q0 and q2 begin to rise (compare e.g. Figs. 3.8a and b). We infer from the results of Chiccoli 
1
et al. [105] that for λ = √ this mismatch is a finite-size effect and that increasing the 
6 1
system size corrects this discrepancy. In our results for λ = √ we choose to interpret the 
6 
temperature at which q0 and q2 are seen to increase as tINB (marked as vertical dashed lines in Fig. 3.8). From Monte Carlo tINB ≈ 1.22, while from mean-field tINB ≈ 1.6. 
q 0 , q 2 
c 

-0.6 

-0.5 
1
-0.8 
-1 
-1 
0 
t 
c ab 
5
1 
0.4 
0.2 
4 
0.5 
0 
3
-0.2 

w

0 

Figure 3.7: Results for λ =0.5 and h =0. Monte Carlo transition temperatures are marked with vertical dashed lines, while the phases are labeled above the plot. a) Monte Carlo (points) and mean-field (lines) uniaxial and biaxial order parameters q0 and q2 (3.21). Note how the change in mean-field q0 is abrupt, indicating the first-order character of the I − NU+ transition. b) Monte Carlo (points) and mean-field (lines) invariant parameter w (1.21) and Monte Carlo energy fluctuation c (3.42) (in green). When calculating w from Monte Carlo simulations results for which Tr(Q2) < 10−3 are discarded, which happens in the isotropic phase. c) Monte Carlo fluctuations of uniaxial (χ0) and biaxial (χ2) order parameters. In Monte Carlo results, at temperatures directly above tINU− ≈ 1.35, minute increase in uniaxiality and biaxiality is visible (see w, q0 and q2), which is identified as a finite-size effect, often observed for such system sizes [104, 105]. The mean-field lines are composed of the high-temperature expansion result for n =8 (see section 3.4.2.2) for t ≥ 1 and of the self-consistent result (see section 3.4.2.1) for t< 1. 
χ 0 q 0 , q 2 
c 

-0.5 

-0.5 
1 -1 
-1 0 
t 
c 
5e-05 
1e-03 4e-05 1e-03 
1e-03 
3e-05 
ab 
5
1
1 
4 
0.5
0.5 
3 

χ2 
w
0 

8e-04 2e-05 6e-04 4e-04 
1e-05 
2e-04 
0 0e+00 
1
Figure 3.8: Results for λ = √ and h =0. Monte Carlo transition temperatures are 
6 
marked with vertical dashed lines, while the phases are labeled above the plot. a) Monte Carlo (points) and mean-field (lines) uniaxial and biaxial order parameters q0 and q2 (3.21). Out of the two equivalent branches of solutions for q0 and q2, the prolate branch (q0 ≥ 0) has been chosen for presentation. b) Monte Carlo (points) and mean-field (lines) invariant parameter w (1.21) and Monte Carlo energy fluctuation c (3.42) (in green). When calculating w from Monte Carlo simulations results for which Tr(Q2) < 10−3 are discarded, which happens in the isotropic phase. c) Monte Carlo fluctuations of uniaxial (χ0) and biaxial (χ2) order parameters. In Monte Carlo results the transition temperature result inferred from the behavior of q0 and q2 differs from what is visible in the fluctuations of energy c and order parameters χ0 and χ2. Peaks in these quantities are much broader than for other phase transitions encountered for the Luckhurst-Romano model (3.1). It can also be told from the behav<BR>ior of order parameter fluctuations χ0 and χ2 (c) that in a broad temperature region below the transition temperature, down to t ≈ 0.7, strong fluctuations are encountered, which affects the precision of determining q0, q2 and w. The mean-field lines are composed of the high-temperature expansion result for n =8 (see section 3.4.2.2) for t ≥ 1 and of the self-consistent result (see section 3.4.2.1) for t< 1. 
λ  tINU (mf)  tINU (MC)  tNU NB (mf)  tNU NB (MC)  tINB (mf)  tINB (MC)  
0.3  1.48  1.08  0.59  0.43  X  X  
0.5  1.85  1.35  0.86  0.62  X  X  
1√ 6  X  X  X  X  1.6  1.22  
Table 3.1: Values of the transition temperatures obtained for the Luckhurst-Romano model (3.1). 



3.6 Studying the model with external field 
We now turn to the general case, described by the model (3.14), which we investigate using the methods described earlier. 
3.6.1 Critical and tricritical points 
The high-temperature expansion mean-field method is useful in finding critical and tricritical points, because it correctly recovers the transition order. A critical or tricritical point can be found by finding a set of conditions for the per-particle free energy (3.33), which are met at such a point exclusively. These conditions can be solved for t = tc and ΔEh2 =ΔEh2 to find the location of the critical or tricritical point in the (t, ΔEh2) plane,
c 
keeping λ as a constant parameter. We find that these conditions attain a particularly compact and simple form if a different decomposition of Q (rather than (3.20)), through the parameters q0 and q2 in the director basis) is adopted. Throughout this section we decompose the alignment tensor Q in a basis fixed on the direction of the field and the laboratory frame of reference {ˆy, ˆ
x, ˆz}: 
(2) (2)
Q = r0T({ˆy, ˆ({ˆy, ˆ(3.46)
x, ˆz})+ r2Tx, ˆz}),
02 
(2)
where zˆh, because we put nˆThe tensors T({ˆy, ˆ=
≡ ˆhz. m x, ˆz}),m 0, 2 are of the form (2.8-2.9), where we substitute {eˆx x, ˆ≡ ˆez z}. The scalar parameters
≡ ˆey y, ˆ≡ ˆr0 and r2 retain the interpretation of the uniaxial and biaxial order parameters, related through linear combinations to q0 and q2, and thus the minimization conditions for the per-particle free energy are similar: 
⎧⎪⎪ 
⎪
⎪
⎨ ⎪
⎪ 
⎪
⎪
⎩ 

∂ ∂r0 fn =  0  
∂ ∂r2 fn = ∂2 ∂r2 0 fn >  0 0  (3.47)  
∂2 ∂r2 2 fn >  0.  

fn is the expanded mean-field per-particle free energy (3.33) with Q expressed using the decomposition (3.46): 
z 
22
(r0 + r2) − t log dΩˆ
n 
−i 1 i 
fn 
= 

t

Ω) , (3.48)
−HMF (ˆ
2 

i! 
i=0 
where: 
(2) (2)
HMF (ˆ0 x, ˆz})+ r2T2 x, ˆz}) · Q(ˆΩ)αβ. (3.49)
Ω) = −E(r0T({ˆy, ˆ({ˆy, ˆΩ) − EΔEhαhβQ(ˆ
3.6.1.1 Mean-field calculation of the critical point 
We proceed to describe the critical temperature and critical field for the critical point of the NP − NU transition upon continuous change of the parameter λ. This is achieved by first observing that at the critical point the parameter r0 has an inflection point. The inflection point at t = tc and h = hc is described by equalities: 
⎧ ⎪
⎪
⎪
⎨ ⎪ 
⎪
⎩ 

∂t  
=  0  
∂r0  
∂2t  
∂r2 0 =  0  (3.50)  
∂3t  
∂r3 0 =  0.  

Furthermore, we assume r2 =0 since at non-zero field the transition is from paranematic to nematic (NP − NU ), and there is no biaxial ordering. Subsequently: 
⎧⎪⎪
⎪
⎨ ⎪
⎪
⎩ 

∂t ∂r0  = ∂t ∂f ∂f ∂r0 = 0  
∂2t ∂r2 0 = ∂ ∂r0  ∂t ∂f ∂f ∂r0  = 0  (3.51)  
∂3t ∂r3 0 = ∂2 ∂r2 0  ∂t ∂f ∂f ∂r0  = 0,  

where f is the free energy per particle. In conclusion, the conditions for the critical point 

can be written as:

⎧ ⎪⎨ ⎪⎩ 

∂f ∂2f∂3f 
=0, =0, =0 
∂r0 ∂r2 ∂r3 
00 (3.52) 
r2 =0. 
Anadditionalrequirementisthat r0 minimizesthefreeenergyperparticle f. Weproceed 
by numerically solving (3.52) for t and ΔEh2, with f ≈ fn 
√ � 
as in the expansion (3.48), in 
3
the order n = 12 for different values of the parameter λ ∈ 0, to obtain the critical
2 
values tc and hc. Because the series in (3.48) involves λ in powers up to n, the series 
16 
converges faster for low λ. Therefore, a more reliable result for λ> 

is obtained by 

first solving (3.52) for λ ∈ [0,

√
1√
] and mapping the outcome to ( 1
66 
, 

3 
2 
] by using the 

duality transformation (3.15). The result is presented in Figs. 3.9a-c. A view of the behavior of q0 upon passing the critical point is presented for λ =0.3 in Fig. 3.10. 
As can be seen in figures 3.9a and 3.9b, the critical field experiences an inflection point <BR><BR>
in the vicinity of λ 

= 

1√ 6
, where hc =0. The presence of the critical point at zero 
2 
c 
ΔEh 

field indicates that the I − NB transition is second-order, which is known to be the case [104, 105]. 
ab 

5.5 
0.02 
5 
t c 

0.01 4.5 
0 
4 
-0.01 
3.5 3
-0.02 
2.5
-0.03 
2 
-0.04 
1.5 
-0.05 
1 
c 

a) ΔEh

2 
c 
versus the biaxiality parameter λ. Note that the critical field decreases when 

molecular biaxiality is increased (i.e. λ 

= 

1√ 6
is approached from either side). Inset: 

1√ 6
(area indicated by a dashed rectangle 

Zoom on the vicinity of the self-dual point λ = 

on the main plot). 

2 
c
b) Critical temperature versus ΔEh

. Points mark the critical point for special values of 

1√ 6
(area indicated 

λ, as labeled. Inset: Zoom on the vicinity of the self-dual point λ = 

by a dashed rectangle on the main plot). 
c) Critical temperature versus the biaxiality parameter λ. 

r 0 

0.3 
0.25 
0.2 
0.15 
0.1 
0.05 0 
ΔEh2 
0.0042 0.0054 0.0066 0.0078 0.0090 0.0102 0.0114 
Figure 3.10: The uniaxial order parameter r0 as defined in (3.46) passing through the critical point at ΔEh2 ≈ 0.007976,tc ≈ 1.4566 for λ =0.3. Note how for larger values of
c 
h the uniaxial order increases in the NP phase above the transition temperature, which is also increased. Above hc, the change in r0 becomes smooth, as the phase transition vanishes. Plots obtained by using the mean-field expansion (3.48) to the order n =6. tc and hc obtained for n = 12. 
3.6.1.2 Mean-field calculation of the tricritical point 
For the NP − NB transition the biaxial order parameter r2 is zero in the NP phase, but is non-zero in the NB phase. Initially discontinuous in the first-order transition, at the tricritical point both r0 and r2 experience an inflection point at temperature tc and field magnitude hc. The corresponding conditions for the mean-field free energy per particle 
at t = tc and h = hc are:
⎧ ⎪
⎪
⎪
⎨ ⎪⎪
⎪
⎪
⎩ 

∂f 
=0 
∂r0 ∂f 
=0 
∂r2 
(3.53)
∂2f ∂r02 =0 
∂2f ∂r22 =0. 
As previously, r0 and r2 need to obey the minimum f condition. We use the expansion 
of the mean-field free energy per particle (3.48) to the order n = 10 and solve (3.53) 
3 
]. As before, we first solve

2
numerically for t and ΔEh2 for different values of λ ∈ [0, 
13
for λ ∈ [0, √ ] and map the result to (√1, ] using the duality transformation (3.15).
2
66 
The result is presented in Figs. 3.11a-c. Exemplary behavior of order parameters as the tricritical point is passed is presented in figures 3.12a-b. 
As we can see in Fig. 3.11a, the tricritical point for λ = 

1√ 6
is present at zero field, at 

tc ≈ 1.6 and coincides with the critical point, thus forming the multicritical point. This further reflects the fact that the I − NB transition is second-order in the self-dual case. 
ab 
c 

t c 
0 
2.5 2 
-0.05 
1.5 

5 
0.1 
4.5 4 
3.5 3 0.05 
c 

a) ΔEh

2 
c 
versus the biaxiality parameter λ. Note that the critical field decreases when 

molecular biaxiality is increased (i.e. λ =

1√ 6
is approached from either side). 

b) Critical temperature versus ΔEh

2 
c 
. Points mark the tricritical point for special values 

of λ, as labeled. 
c) Critical temperature versus the biaxiality parameter λ. 

2 
ΔEh 

r 0 
ΔEh2 
-0.06 
-0.03 
0.2 -0.032 -0.08 
-0.034 -0.036 -0.038 
ΔEh2 0.15 -0.04
-0.1 
-0.042
-0.03 
-0.032 
-0.12 -0.034 
r 2 
-0.036 0.1 -0.048 
-0.038  
-0.14  -0.04  
-0.042 -0.044  0.05  
-0.16  -0.046 -0.048  
-0.18  -0.05  0  

Figure 3.12: Behavior of r0 and r2, as defined in (3.46), when passing through the tricritical point at ΔEh2 ≈−0.0397,tc ≈ 1.4678 for λ =0.3. Plots and values of tc and
c 
hc obtained by using the mean-field expansion (3.48) to the order n = 10. 
a) Uniaxial order parameter r0. Note how for larger values of h the uniaxial order 
increases in the NP phase above the transition temperature, which is also increased. 
b) Biaxial order parameter r2. Note how the biaxial order is not affected by increasing 
h in the NP phase above the transition temperature, which rises with increasing field 
magnitude. 

3.6.1.3 Conclusions and comments on simulations and experiment 
We have arrived at a conclusion that the critical point is not possible to locate in our Monte Carlo setup, because determining the transition order by observing order parameter discontinuities proves infeasible, as the first-order transitions are very weak. Observing the peak in the specific heat (fluctuation of energy) as the field is amplified also gives very little indication as to the position of the critical field, because a) the discontinuity in c due to latent heat at the first-order phase transition is not detectable at the considered temperature resolution, and b) there are no instances of c suddenly becoming very large at some critical field for the considered system size, which would reflect the fact that c is singular at the critical point. Moreover, the peak in c decreases in a smooth manner over a broad range of field magnitudes (see Fig. 3.14 for ΔEh2 > 0). In conclusion, without resorting to e.g. finite-size scaling techniques, for the considered systemsize, MonteCarlo sample sizeandtemperatureresolutionwe areunabletodiscern the order of the transition, and therefore also to locate the critical point itself. 
Inour mean-field resultthe critical point is found forafield magnitudethat is small when compared to the temperature. However, in reality this critical field is very large and difficult to achieve experimentally [113–116], which leads us to expect that the Monte Carlo result should likely reflect the reality of the experiment, i.e. that the critical field should be orders of magnitude higher than mean-field theory predicts. Concerning the tricritical point, as it is predicted to exist for a higher critical field than the critical point, the critical field should be even higher in simulation, although, as is shown later, the structures of the phase diagrams agree for mean-field and simulation results. Interestingly, though it is widely established that both the critical and tricritical points are present for extremely high fields for known mesogens, the behavior of both points predicted by the calculation just presented (Figures 3.9 and 3.11) is such that the critical field decreases as 
1
the self dual point at λ = √ is approached from either side, i.e. as molecular biaxiality 
6 
is increased. Because it appears that the critical field can be made arbitrarily small for molecules which are biaxial enough, this prediction can have experimental implications. 

3.6.2 Results for λ =0.3 
For λ =0.3 we have studied the phase diagram in the (t, ΔEh2) plane using both the mean-field method and Monte Carlo simulations for a range of field magnitudes h ∈ [0, 0.9] for ΔE = ±1. Such field magnitudes can be considered large in relation to e.g. the mean-field result the critical point, which occurs for h2 ≈ 0.008 for λ =0.3. According 
c 
to the discussion conducted in section 3.3.2, for ΔE = +1 we should observe uniaxialitypromoting effects and a critical point for the NP + − NU+ transition, while for ΔE = −1 biaxiality-promoting effects should be encountered, along with a tricritical point for the NP − − NB transition. 
The mean-field phase diagram, displaying first order and second order transition lines, as well as critical and tricritical points, is presented in Fig. 3.13. The mean-field result confirmstheaboveanticipationsandLandau-deGennestheorypredictionscompletelyfor the isotropic-uniaxial nematic transition. In addition, we have found that the transition temperature for the NU+ − NB transition is lowered with increased field magnitude for ΔE = +1, although the change is minute: a ∼ 1% decrease in tNU+NB is noted over an increase from h2 =0 to h2 =0.6. 
FortheMonteCarloresultitismostillustrativetocomparethemean-fieldphasediagram to the energy fluctuation landscape (c, as defined in (3.42)), pictured in Fig. 3.14. The result confirms the mean-field prediction as to the qualitative effects for both cases of ΔE = ±1. As expected, the transition temperatures are lower than in the mean-field result. However, as discussed earlier in 3.6.1.3, the critical and tricritical points cannot be precisely located based neither on the behavior of the energy fluctuation (or specific heat), nor the behavior of order parameters. Instead, we observe the following effects (as labeled in Fig. 3.14): 
• 	
A – The temperature of the NP + − NU+ transition for ΔE = +1 rises as the field is increased and the transition becomes smoother. However, the peak in c does not vanish completely, but is gradually broadened and lowered. The uniaxial order parameter experiences a smooth increase with lowering temperature over the transition for higher field magnitudes (see Figs. 3.15a,c,e). No qualitative change in any of the computed quantities is noted for the mean-field prediction for h2 ≈ 0.008. In conclusion, we are unable to determine whether the critical point 
c 
lies in the range of the considered field magnitudes, or exists for a larger value outside of the studied range. 
• 	
B – The temperature of the NP − − NB transition for ΔE = −1 rises as the field is increased. The peak in c is initially decreased, but stabilizes at c ≈ 4. The NP − − NB transition remains first order and the discontinuity of the uniaxial order parameter q0 is in fact increased when the field is amplified (see Figs. 3.15a,c,e). This behavior is expected to continue until the field is large enough that field-ordering of secondary molecular axes dominates over the ordering of primary molecular axes upon the NP − − NB transition, which will result in a tricritical point. However, for the considered field magnitudes we have not observed this to happen. Furthermore, no qualitative change in any of the computed quantities is noted for the mean-field prediction for h2 ≈ 0.04, which leads us to the conclusion 
c 
that in Monte Carlo simulations the tricritical point is present for a much higher field than the mean-field result would suggest. 
• 	
C – There is no measurable change in the NU+ − NB transition temperature and the peak in c undergoes no systematic changes. We conclude that for the considered field magnitudes this transition is unaffected by increased field in Monte Carlo simulations. 

• 	
D – For non-zero field, the NU+ − NB transition begins to fade away, although this does not occur instantly in the simulation result. The peak in c is gradually smoothed out and moves to higher temperatures with increased field. Eventually, for ΔEh2 ≈−0.25, it is no longer observed. The biaxial order parameter q2 experiences a smooth increase with lowering temperature in place of the zero-field transition (see Figs. 3.15a,c,e). 



Figure 3.13: Mean-field phase diagram of the model (3.14) for λ =0.3. Solid lines represent first-order phase transitions, while second-order transitions are indicated by dashedlines. Criticalandtricriticalpointsareindicatedbyblackfilledcirclesintheinset. The critical point is located at (tc, ΔEh2)cp ≈ (1.4566, 0.007976) and the tricritical point is located at (tc, ΔEh2)tcp ≈ (1.4678, −0.0397) (result obtained by the method explained in section 3.6.1). The high-temperature transition line (NP + − NU+ and NP − − NB) was obtained with the expansion method, described in section 3.4.2.2, with n = 10 in the expression of the approximate per-particle free energy fn (3.33). The NU+ − NB transition line was obtained with the self-consistent method, described in section 3.4.2.1. 

q 0 , q 2 q 0 , q 2 q 0 , q 2 
c 
cc 
0.2 

3 
0 
-0.5 
2 
-0.2 
1 
-0.4 
-1 
0 
tt 
cd 

8
1 
1 

7 
0.8 
6
0.5 
0.6 
5 
0.4 
ab 
8
1 
1 
7 
0.8 
6
0.5 
0.6 
5 
0.4 
0 

www

0 
0.2 
3 
0 
-0.5 
2 
-0.2 
1 
-0.4 
-1 
0 
tt 
ef 

8
1 
1 
7 
0.8 
6
0.5 
0.6 
5 
0.4 

0 
0.2 
3 
0 
-0.5 
2 
-0.2 
1 
-0.4 
-1 
0 
t 
Figure 3.15: Results for λ =0.3 and ΔE = +1. (Caption on adjacent page). 
Figure 3.15: (Left) Monte Carlo results (points) for λ =0.3 and ΔE = +1, compared 
with mean-field results (solid lines), for h =0.3 (a-b), h =0.6 (c-d) and h =0.9 (e-f). 
Monte Carlo transition temperatures are marked with vertical dashed lines, while the 
phases are labeled over the plot. The transition temperatures are estimated based on 
the peaks in the energy fluctuation c. 

Left column (a,c,e): uniaxial and biaxial order parameters q0 and q2, as defined in 
(3.21). <BR><BR>Right column (b,d,f): invariant order parameter w, as defined in (1.21), and energy <BR><BR>fluctuation c, as defined in (3.42). <BR><BR>
Note how with increased field the NP + −NU+ transition shifts to higher temperatures, q0 
increases in the NP + phase, the transition becomes smoother and the peak in c decreases 
and broadens. The NU+ − NB transition appears unaffected by the amplified field. 

The mean-field lines are composed of the high-temperature expansion result for n =8 
(see section 3.4.2.2) for t ≥ 1.4 and of the self-consistent result (see section 3.4.2.1) for <BR><BR>t< 1.4. The dividing temperature has been chosen this high because the mean-field 
free energy per particle expansion (3.33) contains ΔEh2 to the order of n, which impairs 
convergence for h �

=0 in lower temperatures. 
q 0 , q 2 q 0 , q 2 q 0 , q 2 
c 
cc 
0.2 

3 
0 
-0.5 
2 
-0.2 
1 
-0.4 
-1 
0 
tt 
cd 

8
1 
1 

7 
0.8 
6
0.5 
0.6 
5 
0.4 
ab 
8
1 
1 
7 
0.8 
6
0.5 
0.6 
5 
0.4 
0 

www

0 
0.2 
3 
0 
-0.5 
2 
-0.2 
1 
-0.4 
-1 
0 
tt 
ef 

8
1 
1 
7 
0.8 
6
0.5 
0.6 
5 
0.4 

0 
0.2 
3 
0 
-0.5 
2 
-0.2 
1 
-0.4 
-1 
0 
t 
Figure 3.16: Results for λ =0.3 and ΔE = −1. (Caption on adjacent page). 
Figure 3.16: (Left) Monte Carlo results (points) for λ =0.3 and ΔE = −1, compared with mean-field results (solid lines), for h =0.3 (a-b), h =0.6 (c-d) and h =0.9 (e-f). Monte Carlo transition temperatures are marked with vertical dashed lines, while the phases are labeled over the plot. The transition temperatures are estimated based on the peaks in the energy fluctuation c. 
Left column (a,c,e): uniaxial and biaxial order parameters q0 and q2, as defined in 
(3.21). <BR><BR>Right column (b,d,f): invariant order parameter w, as defined in (1.21), and energy <BR><BR>fluctuation c, as defined in (3.42). <BR><BR>
Note how with increased field the NP − − NU− transition shifts to higher temperatures and the peak in c initially decreases, but is retained. The reorientation of the director nˆfrom parallel tonh to perpendicular to ˆh is reflected in the change of sign in q0 and w upon the transition. The NU− − NB transition appears to exist for small values of the field 
''
as a transition between two biaxial phases, N − NB, where N is only distinguished 
BB 
by the fact that its biaxial ordering is field-induced and very small. This transition is shifted towards higher temperatures, but is vanished completely as of h =0.5 both in c and in the behavior of q2. 
The mean-field lines are composed of the high-temperature expansion result for n =8 (see section 3.4.2.2) for t ≥ 1.4 and of the self-consistent result (see section 3.4.2.1) for t< 1.4. The dividing temperature has been chosen this high because the mean-field free energy per particle expansion (3.33) contains ΔEh2 to the order of n, which impairs convergence for h �
=0 in lower temperatures. 

3.6.3 Results for λ =0.5 
1
Because of the self-duality relation (3.15), the effects encountered for λ> √ are equiv
6 1
alent to those present for λ< √ for the opposite cases of the sign of the anisotropy 
6 
of the field polarizability (ΔE →−ΔE) in terms of the structure of the phase diagram and symmetries of the encountered phases. However, since the duality transformation also entails a change of the preferred molecular axis (such that the primary axis is cˆfor 
111
λ< √ and ˆb for λ> √ ), prolate phases (NU+ and NP +) for λ< √ match oblate 
66 6 
1

phases (NU− and NP −) for λ> √ , and vice-versa. This dual symmetry is apparent 
6 
when comparing the mean-field phase diagrams for λ =0.3 (Fig. 3.13) and λ =0.5 (Fig. 3.17). 
According to discussion from section 3.3.2 and Landau-de Gennes theory predictions, for ΔE = +1 biaxiality-promoting effects should be observed, along with a tricritical point of the NP + −NB transition, while for ΔE = −1 one should encounter uniaxiality-promoting effectsandacriticalpointofthe NP −−NB transition. Thisisconfirmedinthemean-field results, as presented in the mean-field phase diagram in Fig. 3.17. Additionally, as in the λ =0.3 case, we have found that the temperature of the NU− − NB phase transition weakly decreases with increased field magnitude, i.e. we observe a ∼ 1% decrease in tNU−NB over an increase from h2 =0 to h2 =0.4. 
ΔEh 2 

0.06
0.3 
0.04 
NB NP + 
0.02
0.2 NB NP +
NP − 
I 
0 
NU− 
0.1 -0.02 
1.82 1.84 1.86 
0 
-0.1 
-0.2 NB NU− 
-0.3 
0.6 0.8 1 1.2 1.4 1.6 1.8 2 
t 

Figure 3.17: Mean-field phase diagram of the model (3.14) for λ =0.5. Solid blue lines represent first-order phase transitions, while second-order transitions are indicated by dashed blue lines. Critical and tricritical points are indicated by black filled circles in the inset. The critical point is located at (tc, ΔEh2)cp ≈ (1.83, −0.004244) and the tricritical point is located at (tc, ΔEh2)tcp ≈ (1.8465, 0.03164) (result obtained by the method explained in section 3.6.1). The high-temperature transition line (NP + − NU+ and NP − − NB) was obtained with the expansion method, described in section 3.4.2.2, with n = 10 in the expression of the approximate per-particle free energy fn (3.33). The NU+ − NB transition line was obtained with the self-consistent method, described in section 3.4.2.1. 
3.6.4 Results for λ = 

6 
For the self-dual case the biaxial nematic phase is of maximum biaxiality, i.e. the alignment tensor Q possesses no purely uniaxial contribution, as its eigenvalues obey the relation ωˆ= −ω ˆ∧ ωˆ=0 and it is proportional to the purely biaxial tensor 
n ml 
1√ 
(2) 
ˆ
m, ˆl)= 
2 
1√
(ˆm ⊗ mˆ− ˆ

n ⊗ nˆ). Any finite finite field introduces a uniaxial contri-

T

({ ˆ

n, 

bution to Q and, depending on ΔE, makes one of the axes cˆand ˆb preferred. Because Q 

is traceless, a field-induced uniaxial contribution results in a decrease in the biaxial part 
of Q. Therefore, the biaxial nematic phase of maximum biaxiality occurs only for h =0. 

b are interchangeable, provided that simultaneously 
6 
we change ΔE →−ΔE (which follows from the duality transformation (3.15)). Such transformation leaves the Boltzmann factor e−H/t for the Hamiltonian (3.14) invariant. 
1√
the molecular axes cˆand ˆ

For λ = 

Therefore, the phase diagram for λ = 

6 
1√
in the (t, ΔEh2) plane is symmetric with respect 
to areflectionaboutthe t-axis. Theisotropicphaseissupersededbyaparanematicphase for finite field magnitudes, NP + for ΔE = +1 and NP − for ΔE = −1. As we can see in the mean-field phase diagram (Fig. 3.18), the I − NB transition temperature drops rapidly upon applying even a small external field. The zero-field phase transition is second-order and the mean-field result suggests that also the finite-field phase transition is second-order as well. However, as of now, phenomenological theory of the NP − NB transition for the self-dual case is not available. 
We present simulation results for small field values of ΔEh2 = ±0.018 in Figs. 3.19-3.19. As we can see, for non-zero field the transition to the biaxial phase is retarded towards lower temperatures (observe q2 in the aforementioned figures and compare with zero-field results in Fig. 3.8), while at t ≈ 1.2 there remains an increase in uniaxial ordering (see q0 in aforementioned figures), which qualitatively is no different from a phase transition in Monte Carlo simulations. Therefore, for the simulation result we distinguish a narrow region where the NU phase is present between NP and NB phases. Comparing the order parameter susceptibilities in Figs. 3.19c and 3.20c, we see that the NP − NU transition is accompanied by an increase in χ2 and the NU − NB transition is accompanied by an increase in χ0. 
As in the result for h =0, there are discrepancies between transition temperatures inferred from increases in q0 and q2 and peaks in c, χ0 and χ2, attributed to a finite-size effect. As in section 3.5, we choose to determine the transition temperatures from the behavior of q0 and q2 for the self-dual case. 

1
√ 
6 
order transition line is indicated by a blue dashed line. The multicritical point is marked by a solid filled black circle at (tc, ΔEh2)mcp ≈ (1.6, 0) (result obtained by the method explained in section 3.6.1). The NP −NB transition line was obtained with the expansion method, described in section 3.4.2.2, with n = 10 in the expression of the approximate per-particle free energy fn (3.33). 
χ 0 q 0 , q 2 
c 

-0.5 

-0.5 
1 -1 
-1 0 
t 
c 
5e-05 
1e-03 4e-05 1e-03 
1e-03 
3e-05 
ab 
5
1
1 
4 
0.5
0.5 
3 

χ2 
w
0 

8e-04 2e-05 6e-04 4e-04 
1e-05 
2e-04 
0 0e+00 
1
Figure 3.19: Results for λ = √ and ΔEh2 = −0.018. Monte Carlo transition tempera
6 
tures are marked with vertical dashed lines, while the phases are labeled above the plot. a) Monte Carlo (points) and mean-field (lines) uniaxial and biaxial order parameters q0 and q2 (3.21). Out of the two equivalent branches of solutions for q0 and q2, the prolate branch (q0 ≥ 0) has been chosen for presentation. b) Monte Carlo (points) and mean-field (lines) invariant parameter w (1.21) and Monte Carlo energy fluctuation c (3.42) (in green). When calculating w from Monte Carlo simulations results for which Tr(Q2) < 10−3 are discarded, which happens in the isotropic phase and in the NP + phase when the field-induced ordering is very small. c) Monte Carlo fluctuations of uniaxial (χ0) and biaxial (χ2) order parameters. 
Compare the presented results with Fig. 3.8 to note that the transition to the NB phase shifts considerably towards lower temperatures (tNP −NB ≈ 0.98) when even a small field is applied. In Monte Carlo results we also distinguish an intermediate nematic uniaxial (NU−) phase, characterized by substantially larger uniaxial ordering than in the fieldinduced paranematic phase (NP −), where it is below measurement precision (for the considered field magnitude). 
The mean-field lines are composed of the high-temperature expansion result for n =8 (see section 3.4.2.2) for t ≥ 1 and of the self-consistent result (see section 3.4.2.1) for t< 1. 
q 0 , q 2 
c 

ab 
5
1
1 
4 
0.5
0.5 
3 
-0.5 

-0.5 
1 
-1 
-1 
0 
t 
c 

w

0 

1
6 
are marked with vertical dashed lines, while the phases are labeled above the plot. a) Monte Carlo (points) and mean-field (lines) uniaxial and biaxial order parameters q0 and q2 (3.21). Out of the two equivalent branches of solutions for q0 and q2, the prolate branch (q0 ≥ 0) has been chosen for presentation. b) Monte Carlo (points) and mean-field (lines) invariant parameter w (1.21) and Monte Carlo energy fluctuation c (3.42) (in green). When calculating w from Monte Carlo simulations results for which Tr(Q2) < 10−3 are discarded, which happens in the isotropic phase and in the NP + phase when the field-induced ordering is very small. c) Monte Carlo fluctuations of uniaxial (χ0) and biaxial (χ2) order parameters. 
Compare the presented results with Fig. 3.8 to note that the transition to the NB phase shifts considerably towards lower temperatures (tNP +NB ≈ 0.98) when even a small field is applied. In Monte Carlo results we also distinguish an intermediate nematic uniaxial (NU+) phase, characterized by substantially larger uniaxial ordering than in the fieldinduced paranematic phase (NP +), where it is below measurement precision (for the considered field magnitude). 
The mean-field lines are composed of the high-temperature expansion result for n =8 (see section 3.4.2.2) for t ≥ 1 and of the self-consistent result (see section 3.4.2.1) for t< 1. 


3.7 Summary 
We have explored the effects of a general vector field in the l =2 dispersion model of biaxial liquid crystals (3.14). Three distinct classes of phase diagrams in the (t, ΔEh2)plane have been identified for the prolate, oblate and self-dual case. Examples have been studied in detail with mean-field methods for the prolate case of λ =0.3 (Fig. 3.13), for the oblate case of λ =0.5 (Fig. 3.17) and for the self-dual case (Fig. 3.18). The behavior of the critical and tricritical points has been studied for the entire spectrum of 
3
the biaxiality parameter λ ∈ [0, ], presented in Fig. 3.9 and Fig. 3.11. 
2 
3.7.1 Phase diagrams and order parameters 
3.7.1.1 λ =0.3 
For the prolate case of λ =0.3, mean-field calculations show an increase in the I − NU+ phase transition temperature with increased field, while the isotropic phase I is replaced by the field-induced paranematic prolate phase NP + for positive sign of the anisotropy of the field polarizability ΔE = +1 and by the paranematic oblate phase NP − for ΔE = −1. The NP + − NU+ phase transition, occurring for ΔE = +1 is first-order, until it is terminated at a critical point for critical field hc ≈ 0.008. For ΔE = −1, field-induced biaxiality emerges in the zero-field nematic uniaxial prolate phase NU+. The NP − − NB transition is initiallyfirst-order, but changestosecond-orderasthecriticalfield hc ≈ 0.04 is passed. The NU+ − NB transition, occurring for ΔE = +1 is largely unaffected by the increased field, however we do note a ∼ 1% decrease in the transition temperature, when the field is increased from h2 =0 to h2 =0.4. For ΔE = −1 said transition disappears, as both phases are now biaxial. 
By Monte Carlo simulation we have obtained the energy fluctuation landscape as a function of t and ΔEh2 for the case of λ =0.3 for h ∈ [0, 0.9] and ΔE = ±1, presented in Fig. 3.14. The result shows qualitative agreement with mean-field predictions, presented in Fig. 3.13, as to the structure of the phase diagram and the phases encountered. However, the methods employed did not allow for precisely determining the position of the critical point, which appears to lie for a field magnitude much higher than mean-field theory predicts. The tricritical point is believed to exist for a critical field well outside the considered range of field magnitudes. 
Order parameters from both Monte Carlo simulations and mean-field calculations have been compared in Figs. 3.15-3.16. The discussed effects are identified by the behavior of order parameters q0 and q2 (3.21) and the invariant order parameter w (1.21). Namely, in paranematic phases the uniaxial order parameter q0 is non-zero and increases as the field is amplified, while its sign identifies whether the phase is of the prolate or oblate variety (q0 > 0 for NP + and q0 < 0 for NP −). For ΔE = +1 the gradual smoothing of the NP + − NU+ phase transition is evidenced by the smoothing of q0 and the peak in the energy fluctuation c, defined in (3.42). For ΔE = −1, q0 undergoes a change of sign from q0 < 0 to q0 > 0 upon the NP − − NB transition, which reflects the fact that in the NP − phase the field orients the director nˆparallel to the field, while upon the phase transition nˆis aligned in perpendicular to the field. As the field is amplified, the biaxial ordering in the field-induced NB phase rises, identified by rising (moduluswise) q2. The zero-field NU+ − NB phase transition vanishes instantly in the mean-field result for h �
=0, however in Monte Carlo results there briefly exists a transition between 
''
two biaxial phases: N − NB, where N is only distinguished by the fact that it is 
BB 
field-induced and has very small biaxial ordering. As the field is increased, this phase transition moves towards higher temperatures and is vanished as of h =0.5. 
3.7.1.2 λ =0.5 
Fortheoblatecasethedualsymmetryofthemodel(3.14),definedby(3.15),ismanifested when one compares the mean-field phase diagram for λ =0.5 (Fig. 3.17) with the phase diagram for the prolate case for λ =0.3 (Fig. 3.13). Namely,thephasediagramsareones mirror images (qualitatively), with the NP + and NU+ phases exchanged with NP − and NU− phases respectively. To avoid redundancy, the discussion of the effects is omitted. 
3.7.1.3 λ = √1
6 1
The self-dual case, λ = √ , for h =0 contains a biaxial phase which is of maximum 
6 
biaxiality, i.e. the alignment tensor Q does not contain a purely uniaxial contribution. Field-induced uniaxiality for h �
=0 causes a decrease in biaxiality and lowers the I − NB transition temperature rapidly, even for small fields. The high temperature phase has field-induced uniaxial ordering and is therefore a paranematic phase, NP + and NP − for ΔE = +1 and ΔE = −1 respectively. We have obtained a mean-field phase diagram for 
1
λ = √ , presented in Fig. 3.18. Because of dual symmetry, defined in (3.15), the phase 
6 
diagram is symmetric with respect to a reflection about the t-axis. There is only one transition line in the diagram, the NP − NB phase transition line, with a multicritical point at (tc,hc) ≈ (1.6, 0), where six phases: I,NU+,NU−,NP +,NP −,NB meet (when considering both (t, ΔEh2) and (λ, t) phase diagrams). 
1
The effects of an external field on the behavior of order parameters in the case of λ = √ 
6 
arepresentedinFig. 3.19-3.20,wherebothmean-fieldandMonteCarlosimulationresults are presented. In both cases of ΔE = ±0.018 the Monte Carlo transition to the NB phase is retarded to tNP NB ≈ 0.98 with respect to the zero-field transition at tINB ≈ 1.22, as presented in Fig. 3.8. In Monte Carlo results we are compelled to treat the increase in q0 at t ≈ 1.2 as a separate phase transition, therefore there is a region where in simulations an intermediate uniaxial nematic phase exists, NU+ for ΔE = +1 and NU− for ΔE = −1, between NP and NB phases. 
3.7.1.4 Mean-field study of critical and tricritical points 
The mean-field high-temperature expansion method has been used for a detailed study 
3
of the behavior of critical and tricritical points, depending on the parameter λ ∈ [0, ].
2 
To this end conditions for the critical (3.52) and tricritical (3.53) points have been found and solved numerically for the critical and tricritical field and temperature, with the result presented in Fig. 3.9 and Fig. 3.11. We have taken advantage of the duality transformation (3.15) to improve the convergence of the mean-field expansion of the per<BR>
particle free energy (3.48) by transforming the result for λ ∈ [0,

√
1√
] to λ ∈ ( 1
66 
, 

3 
2 
]. 

The result confirms that the critical and tricritical points are located for ΔEh2 > 0 for 
√
1166 
√
and ΔEh2 respectively and travel towards ΔEh2 =0 as λ = 
1√ 
6 
< 0 for λ> 

λ< 

is approached from either side, that is, the critical field is seen to decrease as molecular 

biaxiality is increased. At λ 

= 

1√ 
6 
both points coincide, thus forming the multicritical 

point, at which, in the “super” phase-diagram in the parameter space of (t, λ, ΔEh2), six 
phases meet: I,NU+,NU−,NP +,NP − and NB. As λ is increased over 
1√
, the critical 
6 
and tricritical points emerge at opposite sides of the ΔEh2 =0 line, for ΔEh2 < 0 and ΔEh2 > 0 respectively, which reflects the dual symmetry of the model (3.14). 
Chapter 4 



Tetrahedratic order and chiral symmetry breaking in a dispersion model of liquid crystals 
4.1 Introduction 
Developments of the last two decades show that bent-shaped mesogens (bent-core compounds, certain ferrocene mesogens and flexible dimers) exhibit many fascinating structures [36, 37, 41, 44, 117]. Of particular interest is the possible existence of biaxial order [27, 29, 30, 32, 88], spontaneous chirality and formation of chiral supramolecu<BR>lar structures and helices [38–40, 43], the twist-bend nematic phase [45–50], as well as tetrahedratic order [51, 52, 54, 55]. In the present chapter we explore a minimal-coupling dispersion model [1, 118] which exhibits biaxiality, tetrahedratic order and chiral symme<BR>try breaking and can be used to describe behavior in certain systems of supramolecular bent mesogen complexes. The results presented in this chapter have been published in [119]. 
The minimal-coupling model is devised upon the arguments of Lubensky and Radzihovsky [51], who give thorough group-theoretical classification of tensorial order param<BR>eters and phase sequences which can emerge in bent-core systems. Their finding shows that to account for chiral ordering, at the microscopic level interactions should involve couplings of biaxial quadrupolar (second-rank) and tetrahedratic octupolar (third-rank) tensors, as well as a vector to account for the static (steric or electric) dipole. The inclusion of tetrahedratic tensors is also shown to lead to e.g. ambidextrous chirality [55, 120], i.e. twisted domains of opposite handedness. A study involving a purely tetra<BR>hedratic interaction was notably performed by Fel [121], predicting a second-order phase 
81 

transition from isotropic to a nematic tetrahedratic phase. A Monte Carlo study was performed by Romano [122], in which case a weakly first-order transition was found. 

4.2 	Dispersion model with tetrahedratic coupling (infinite pitch limit) 
The model proposed in [1] is designed considering biaxial and tetrahedratic multipolar couplings in the pair dispersion potential: 
(3) (3)
V (Ωˆi, Ωˆj)= −E Q(Ωˆi) · Q(Ωˆj)+ τT(Ωˆi) · T(Ωˆj) . (4.1)
22 
The full Hamiltonian is obtained by performing a sum of terms (4.1) over nearest neigh<BR>bors on a simple cubic lattice of N particles with periodic boundary conditions: 
H  =  1 2  N  V (ˆΩi, ˆΩj).  (4.2)  
(i,j)  

This model corresponds to the general model (2.11) for the case when all parameters aside from λ, E and τ are zero. In our discussion we will assume that the model is studied with respect to the dimensionless temperature t = kBT/E. 
As in the previous cases, the dot products in (4.1) are understood as full contractions over Cartesian indices. The quadrupolar tensor Q(Ωˆi) is defined in (3.2). We recognize that when τ =0, the model reduces to the l =2 model of biaxial nematics (3.1), as studied previously [64, 104, 105] and in Chapter 3, while the term proportional to Eτ is related to the model due to Fel [121]. While the two terms by themselves are rotationally invariant, a potential like (4.1)needsto establish arelationbetweentheorientationofthe 
(2) (3)
axes of Tm and T2 . In (4.1), the quadrupolar part remains identical to (3.1), while in the tetrahedratic term the tensor T(3)(Ωˆi) is the l =3 irreducible spherical tensor (2.10), 
2 
symmetrized with respect to the full tetrahedratic group Td with the assumption that the tetrahedratic twofold symmetry axes coincide with the D2h twofold symmetry axes of the quadrupolar tensor. The potential (4.1) does not contain cross-coupling terms between quadrupolar and octupolar tensors, but it can be shown that this precise arrangement of axes realizes the minimum of the lowest-order cross-coupling terms, which turns out to be zero. 
The choice of symmetry axes has a physical interpretation. The relative alignment of the quadrupolar and tetrahedratic tensors corresponds to a ground state configuration of two bent molecules, in which they align in an embrace-locked fashion, as presented 

ˆ
pair complexes of opposite chirality. The molecular bases {ˆˆc} and {aˆm,bm,cˆm}
a, b, ˆcorrespond to the ordinary and reflected bases of opposite handedness. (Reprinted with permission from [1], �c2009 by the American Physical Society.) 
in Fig. 4.1. As we see, the complexes are non-polar, because the dipoles (parallel to the bisector of the opening angle) cancel out, which justifies the omission of a vector in (4.1), and exhibit conformational chirality as the relative twist angle δ ∈ [0,π) is varied. 
π
For any δ �0 and δ = , the conformation has D2 non-polar chiral symmetry and is 
2 
available in two states, each obtained by a reflection about the origin of the other (left 
π
and right portions of Fig. 4.1). The special cases of δ =0 and δ = correspond to 
2 
non-chiral biaxial D2h and full tetrahedratic Td symmetry, respectively. It is speculated that forming of such complexes as in Fig. 4.1 is responsible for homogeneous chirality and tetrahedratic order [1]. 
Although the discussion just cited is valid for the type of construction pictured in Fig. 4.1, the two states of opposite chirality, parametrized by a binary microscopic variable pi (2.15), may be more general. For instance, flexible dimers and bent-core molecules are terminated with flexible alkyl chains at both ends and have an unlimited number of conformations, most of which are chiral and on average can be divided into two classes of conformations of opposite chirality, parametrized by pi. To both of these classes there may correspond many conformational states of the molecules, but we assume that apart from their chirality they do not modify the quadrupolar interaction in an essential way. Therefore the considered model (4.1) is not limited to a picture like in Fig. 4.1. 
ˆ
In (4.1), chirality manifests itself when reflections of the molecular basis Ωˆ{ˆb, cˆ}
= a, that do not preserve handedness are considered. Under such reflection the quadrupolar tensor Q(ˆ(3) Ω) changes sign. Although 
Ω) is left invariant, while the octupolar tensor T(ˆ
2 a reflection followed by a 90◦ rotation about one of the twofold axes leaves T(3) Ω) 
(ˆ
2 
invariant, it is not a symmetry of Q(ˆThus, the two chiral states (Fig. 4.2) are 
Ω). distinguished by the parity of the molecular basis p = aˆ· (ˆb × cˆ) = ±1, which can be separated as an independent, discrete degree of freedom. Therefore, from now on we rename {ˆˆb, ˆΩ ' and imply that the generalized coordinates for molecule i are Ωˆi ≡
a, c}≡ ˆ

Ω) is visualized in blue; the octupolar tensor Tˆ
2 
visualized in orange. Surfaces obtained as r = exp(Q(ˆu ⊗ ˆ2 Ω) ·
Ω) · ˆu) and r = exp(T(3)(ˆuˆ⊗ uˆ⊗ uˆ), where uˆis a unit vector expressed in polar coordinates. 
ˆ
{pi, (ˆai,bi,cˆi)}. Being an independent degree of freedom, parity needs to be included in the trace in the definition of the canonical ensemble average (1.3), therefore it is assumed that Tr ≡ pi=±1 dΩˆ' i. Taking this degree of freedom into account means that we 
{Ωˆi,pi}
allow the molecules to have two different conformations of opposite chirality. 
A parity flip amounts to inverting one or more of the molecular axes or, equivalently, 
(3)(ˆ
changing the sign of TΩ) (while Q(ˆΩ) remains invariant). Therefore, we can write 
2 
the tetrahedratic term of (4.1): 
(3) (3) (3) (3)
VT (Ωˆi, Ωˆj)= −EτT(Ωˆi) · T(Ωˆj)= −EτT({pi, (ˆai, ˆbi,cˆi)}) · T({pj, (ˆaj, ˆbj,cˆj)})= 
222 2 
(3) (3)
ˆˆ
= −EτpipjT2 ({aˆi,bi,cˆi}) · T({aˆj,bj,cˆj}).
2 
(4.3) Here pi,pj = ±1 areparitiesatthelatticesites i and j. Weobservethatthetetrahedratic term is equal to the scalar product of tetrahedratic tensors in the canonical (righthanded) molecular bases, multiplied by the parities pipj. The minimum of (4.1) with respect to the separated degrees of freedom is such that not only all molecular axes of molecules i and j are in parallel, but also that pi = pj. Thus, the ground state favors homogeneous chirality. However, as we discuss in Chapter 5, the considered symmetries 
(3) (2) (2)
allow for constructing an antisymmetric tensor out of T(ˆ(ˆ(ˆ
Ω), TΩ) and TΩ), and 
20 2 
therefore the ground state is unstable with respect to spontaneous formation of helical twists. Thus, the discussion of this chapter is to be treated as the infinite pitch limit of the proper minimal coupling model. 
4.2.1 Order parameters 
The model (4.1) leads to a multitude of phases of non-standard symmetries. As in the case of the l =2 model of biaxial nematics (3.1), which is responsible for uniaxial and biaxial nematic ordering, the appropriate order parameters for monitoring D∞h and D2h symmetries are the uniaxial and biaxial order parameters q0 and q2, as defined in (3.21), while the definition of the alignment tensor Q (3.17) holds. 
Monitoring of tetrahedratic order is performed through the ensemble average of the Td symmetry-adapted function: 
(3) (3) (3)
Δ= Tm, ˆΩ), (4.4)
({ˆl, ˆn}) · T(ˆ
222 2 
where we have substituted Ωˆ1 ≡{ˆl, m,ˆnˆ} and Ωˆ2 ≡ Ωˆin the definition (2.6) and applied 
(3)
the ensemble average to both sides, assuming that T({ˆl, ˆn}) is constant. Further
2 m, ˆmore, calculating Δ(3) can be done without determining the director basis at all (thus 
22 
limiting the error due to diagonalization), by observing that: 
(3) (3)
(ˆ(ˆ
T2 Ω) · T2 Ω) = 
2 
(3) Ω) · T(3) n}) ⊗ T(3) (3) (3) (4.5) 
= T({ˆm, ˆl, ˆn}) · TΩ) = ,
(ˆl, ˆ({ˆm, ˆ(ˆΔ
2 mm 2 22 
m=0,±1,±2,±3 
where we have inserted the Td-symmetrized identity operator for the l =3 subspace of angular momentum (2.7) and observed the orthogonality relations for irreducible tensors. Thus, we use the identity: 
(3) (3) (3)
(ˆ(ˆ
Δ= TΩ) · TΩ) (4.6)
222 2 
(3)
to calculate Δ
22 . 
Phase chirality is measured by the ensemble average of the molecular parity: 
p = aˆi · (ˆbi × cˆi). (4.7) 
It is shown that the leading mean-field contribution to p reads [1]: 
√ 
2 (3) 3 
p ≈ τλ(−3+2λ2)Δ22 q2, (4.8)
210t4 
which indicates that macroscopic chirality is induced by non-zero tetrahedratic order coupled with biaxiality. (p =0 in uniaxial states.) 

4.2.2 Phases of model (4.1) 
Apart from the phases known to exist for the quadrupolar model (3.1), which are the uniaxial nematic NU and biaxial nematic NB phases, there are additional phases in (4.1), 
in which tetrahedratic order is non-zero: 
• 	The tetrahedratic phase T , characterized by global symmetry of the full tetrahedratic group Td (Fig. 4.3a). This phase is identified by non-zero Δ(3) and vanishing 
22 
uniaxial and biaxial order parameters. The molecular arrangement which produces this phase is achieved when all of the molecular tetrahedratic threefold symmetry axes align in parallel. The relative arrangements of molecular bases which satisfy this condition are realized randomly and on average none of the twofold axes of the quadrupolar tensors are distinguished. 
• 	The nematic tetrahedratic phase NT , characterized by global dihedral symmetry described by the point group D2d (Fig. 4.3b). An example of an object which bears the D2d symmetry is a tetrahedron stretched or compressed along one of its twofold axes. NT differs from T in the respect that also the uniaxial order parameters q0 and Δ(2) are non-zero. This phase is also referred to as “uniaxially deformed 
02 
tetrahedratic phase”. Much like NU , NT comes in two varieties, depending on the type of the deformation: prolate NT + (stretched) and oblate NT − (compressed). The molecular arrangement which produces this phase is achieved when all of the molecular tetrahedratic threefold symmetry axes align in parallel, with an addition that the uniaxial symmetry axes of the quadrupolar part are also parallel. The relative arrangements of molecular bases which satisfy this condition are realized randomly and none of the secondary twofold axes of the quadrupolar tensors are distinguished. 
• 	The chiral nematic tetrahedratic phase N∗ , characterized by global dihedral sym-
T 
metry of the D2 point group (Fig. 4.3c). An example of an object which bears this symmetry is a tetrahedron which is twisted along one of its twofold axes. In this phase, apart from non-zero tetrahedratic and uniaxial order parameters, also 
(2)
the biaxial order parameters q2 and Δ20 are non-zero and, consequently, also is the macroscopic parity p. This is the lowest-symmetry phase that is available in the model (4.1). The molecular arrangement which produces this phase is achieved when all of the corresponding molecular basis vectors are in parallel, i.e. aˆi aˆj, ˆbi ˆbj, cˆi cˆj, while pi = pj. 

All the phase transitions, except for those where the uniaxial parameter assumes a non-zero value, i.e. I − NT (isotropic to nematic tetrahedratic), T − NT (tetrahedratic to nematic tetrahedratic) and I − NU (isotropic to nematic uniaxial), are second-order. The phase diagram in the (λ, t) plane depends greatly upon the parameter τ, which has several special values, with respect to which four classes are distinguished [1, 118, 119]. It has been found via bifurcation analysis that there exists a highly-symmetric value 
28 1
of τ = , for which at the biaxial self-dual point λ = √ and temperature t =8/5 a
15 6 
multicritical Landau point is predicted, where six phases, I, NU−, T , NT +, NT − and NT ∗ , meet. There are two phase diagram classes for τ< 2815 and two for τ> 28 . For τ � 1.54,
15 
the high-temperature regime is dominated by nematic phases (NU , NT and NB) and there is no purely tetrahedratic phase T . For 1.54 � τ< 28 the T phase is found and 
15 
a first-order T − NT + transition is present. For 28 >τ � 6.16 there is no stable NB
15 
phase, while the remaining nematic uniaxial phase is NU−; the I − NU− transition line meets with the NU− − NT −, I − T and T − NT − transition lines at a bicritical point. The last class of diagrams for τ � 6.16 exhibits a direct I − T second-order transition and no purely nematic phases. In all phase diagram classes the temperature of the N∗ 
T 
1
phase is highest at the biaxial self-dual point λ = √ . 
6 

4.2.3 Scope of presented research 
The results presented in subsequent sections pertain to the nematic-dominated phase diagram class for τ � 1.54, i.e. τ =1 in our case, and the phase diagram containing the 
28
multicritical point for τ = , presented in Fig. 4.4b. For these two values Monte Carlo 
15 
temperature scans are performed for two generic values of the biaxiality parameter λ: 
1
0.3 and 0.5, and the biaxial self-dual value λ = √ . In the case of τ =1, the results 
6 
supplement the phase diagram (Fig. 4.4a), obtained in [1], with temperature scans. Throughout the following sections we assume E =1. 

28
15 
1
presence of the biaxial phase in the vicinity of λ = √ . Observe that in b) the T phase
6 
exists in a very narrow temperature and that there exists a direct I − N∗ transition for 
T 
1
λ = √ . 
6 1
(a obtained by Longa et al. [1]. b – points for λ =0.3, 0.5, √ obtained by Monte Carlo 
6 
simulation as presented in section 4.5, remaining points and mean-field lines obtained in collaboration [119]. (a) reprinted with permission from [1], �c2009 by the American Physical Society. (b) c
�2012 by the American Physical Society) 


4.3 Monte Carlo simulations 
In the presented results a simple cubic lattice of N = 16 × 16 × 16 sites with nearestneighbor interactions and periodic boundary conditions is considered. The temperature scan progresses by heating, i.e. the initial lattice is prepared in a completely-ordered N∗ 
T 
stateinthelowesttemperatureandthefinalstateispasseddownasaninitialstateforthe subsequent simulation for a larger temperature. At the beginning of the thermalization processtheradiusoftherandomwalkinrotationspace r (asdefinedin(1.49)iscomputed to match an acceptance rate in the 40% ∼ 50% range. The thermalization consists of 500000 lattice sweeps. Production (equilibrium state sampling) consists of 500000 lattice sweeps, while every tenth configuration is considered for calculating averages to minimize correlations between subsequent states. After the final state is passed onto the next simulation for a higher temperature, another round of 500000 thermalization cycles follows, etc. The ensemble averages and susceptibilities are calculated in the same manner as it was introduced in Chapter 3, specified in the definitions (3.38) and (3.39). 
Alatticesweep consists ofperformingtrialmoveson N randomly chosenmolecules, while a Monte Carlo trial move consists of separately applying the Metropolis algorithm (see section 1.8.1) to the orientational and parity degrees of freedom. Namely, considering a trial orientational move {aˆi, ˆbi,cˆi}→{aˆi' , ˆbi' ,cˆ' } for a particle i, the energy difference for 
i
the potential (4.1) reads: 
z=6
1 ˆˆ' ˆb '' ˆ
ΔEi = V ({pi, (ˆai,bi,cˆi)}, {pj, (ˆaj,bj,cˆj)}) − V ({pi, (ˆai, i,cˆi)}, {pj, (ˆaj,bj,cˆj)}) . 
2 
i(j) 
(4.9) The energy difference upon a parity flip pi →−pi reads: 
z=6 
1 ˆˆˆˆ
ΔEi = V ({pi, (ˆai,bi,cˆi)}, {pj, (ˆaj,bj,cˆj)}) − V ({−pi, (ˆai,bi,cˆi)}, {pj, (ˆaj,bj,cˆj)}) . 
2 
i(j) 
(4.10) Note that there is no adjustable parameter which allows for an adjustment of the acceptance rate for parity moves. We could, in principle, limit the number of accepted moves by reducing the number of attempts, which is not necessary. However, we have no means of encouraging more parity moves to be accepted, which results in parity being frozen for long periods of time in low temperatures, as the relative energy difference associated with a single flip is becoming large. Because the starting configuration is homogeneously chiral, this problem does not affect the lowest temperatures, as the initial state is already close to equilibrium. In intermediate temperatures, however, the time scales for the dynamics of parity are much larger than for rotational degrees of freedom, which causes the Monte Carlo sample for quantities linked to parity to experience non-negligible temporal correlations. 
4.3.1 Susceptibilities 
We investigate the previously defined uniaxial and biaxial order parameter susceptibilities (or measures of order parameter fluctuations) (3.44) and (3.45), but also consider susceptibilities of tetrahedratic order and parity. They are defined analogously, in terms of instantaneous lattice averages, as in the definition (3.39): 
22
χp = p˜M − p˜M , (4.11) 
˜(3) ˜(3) ˜(3) ˜(3) 2
χ (3) = T· TM − T· TM . (4.12)
22 22
T
2 

4.3.2 Correlation functions 
Becauseofthesuspicionthatinlowertemperaturestheremightexistdomainsofopposite homogeneous chirality, the average parity (4.7) can also be obtained by considering a correlation function: 
|i−j|→∞
2
Gpp(|i − j|)= pipj −→ pi · pj = p. (4.13) 
Because of the periodic boundary conditions, the limit for the argument of (4.13) is half of the edge of the lattice L. By calculating p as: 
p = Gpp(L/2), (4.14) 
we evade the error due to spatial correlations when calculating the lattice average parity. 

4.3.3 Determining transition temperatures 
Themostpronouncedphasetransitionsof(4.1)aretheonesassociatedwithtetrahedratic ordering, i.e. the I − T (isotropic to tetrahedratic, not observed in simulations), I − NT (isotropic to nematic tetrahedratic, see Figs. 4.8 and 4.10), I − N∗ (isotropic to chiral 
T nematictetrahedratic,seeFig. 4.9), NU −NT (nematicuniaxialtonematictetrahedratic, see Figs. 4.5 and 4.7) and NB − N∗ (nematic biaxial to chiral nematic tetrahedratic, 
T 
see Fig. 4.6) transitions, which produce a very strong peak in the energy fluctuation c (defined in (3.42)). Recall that the fluctuation of energy calculated for the biaxial model (3.1) was of the order 5 ∼ 7 at the I − NU phase transition for a similar system (see e.g. Fig. 3.6b). The same method of calculating c yields peaks of the order of 30 for the aforementioned transitions in (4.1) (see e.g. Fig. 4.6d). This often means that the peak at the I − NU and I − NB transitions is overshadowed by the peak at a tetrahedratic symmetry breaking transition, if they are close. Therefore, for these transitions we determine the temperature by peaks in χ0 and χ2 respectively, adjusting to agree with the behavior of order parameters. 
The transitions involving spontaneous breaking of tetrahedratic symmetry produce pronounced peaks in c and χ (3) , which provide excellent indication of the transition tem-
T
2 
perature. 
The most challenging are the transitions where spontaneous breaking of chiral symmetry occurs separately from breaking of tetrahedratic symmetry, i.e. NT − N∗ (T phase is 
T 
not observed in our results), observed for τ =1,λ =0.3 (Fig. 4.5), τ =1,λ =0.5 (Fig. 
28 28
4.7), τ = ,λ =0.3 (Fig. 4.8) and τ = ,λ =0.5 (Fig. 4.10). Firstly, we fail to 
15 15 
detect any peak in c due to the transition. Secondly, the collection of equilibrium states from which the averages and susceptibilities are calculated is significantly correlated with respect to parity due to large time-scales, resulting in lowered precision of determining p and χp. If a peak in χp is present and its position is unambiguous, we use it to determine the transition temperature. There are, however, cases where χp does not exhibit an unequivocal peak and in these cases we resort to approximating the transition temperature as the point, where p starts to increase in the N∗ phase. Furthermore, q2
T 
is linked to p through (4.8), which affects the behavior of the order parameters and χ2 in NT and N∗ phases as well. 
T 


4.4 Results for τ =1 
This result is complementary to the phase diagram in Fig. 4.4a. The original result is supplemented in this section by temperature scans presenting thermal behavior of order parameters q0, q2, p and Δ(3), energy fluctuation c and susceptibilities χ0, χ2, χp
22 
and χ (3) . The results are presented in figures 4.5-4.7, where Monte Carlo transition 
T
2 
temperatures are indicated. Below will follow a general description of the results, while detailed comments on the behavior of the measured quantities are given in the captions 
1
of figures: Fig. 4.5 for λ =0.3, Fig. 4.6 for λ = √ and Fig. 4.7 for λ =0.5. 
6 
1
The only instance of NB phase exists for λ = √ (Fig. 4.6), where we have a direct 
6 second-order I − NB transition, evidenced by and increase in q0 and q2, while p and Δ(3) 
22 parameters are zero. A peak in χ2 is also noted and used to determine the transition temperature tI−NB ≈ 1.2. A second-order phase transition to the N∗ phase is subse-
T 
quently encountered as the temperature is lowered, evidenced by an increase in Δ(3) and
22 
peaks in χ (3) and χp. 
T
2 
For the generic prolate (λ =0.3, Fig. 4.5) and generic oblate (λ =0.5 4.7) cases, as the temperature is lowered, first a uniaxial nematic phase occurs, followed by a second-order transition to the nematic tetrahedratic NT phase. Because of the slow dynamics of parity and the fact that biaxial ordering in the presence of tetrahedratic order is linked to parity, p and q2 experience fluctuations due to long temporal correlations in the NT phase, evidenced by the behavior of χp and χ2. The correlations continue to the N∗ 
T 
phase, which is visible in χp, particularly in the case of λ =0.3, however χp indicates that these effects are less pronounced in low temperatures. The transition to the N∗ 
T 
phase is second-order and is evidenced by an increase in p, q2 and a peak in χp, which was used to determine the transition temperature. 
ab 

χ 0 
q 0 , q 2 
χ p χ (3)
T 
2 
1 

1 

0.0014 
0.8 0.0012 
0.8 
0.6 
0.001 
0.60.4 0.0008 
)(3
Δ 
22 , p 
0.2 
0.0006
0.4 
0 
0.0004 
0.2
-0.2 
0.0002 
-0.4 
0 
0 
cd 

7e-05 

0.00025 25 

0.004 
0.0035
6e-05 
0.0002 20 
0.003
5e-05 
c

0.0025
0.00015 15
4e-05 
0.002 
3e-05 
0.0001 
10 
0.0015 
χ2 
2e-05 
0.001 
5e-05 
5 
1e-05 0.0005 
0 0 0
0 
Figure 4.5: Monte Carlo results for τ =1 and λ =0.3. Obtained transition temperatures: tI−NU+ ≈ 1.18,tNU+−NT + ≈ 1.18,tNT +−N∗ ≈ 0.76 (marked as dashed vertical 

T 
lines). Note that the NU+ phase exists in a very narrow temperature range, hence the 
peak in c (d) is not discernible from the NU+ − NT + peak. Peaks of χ0 (c) at I − NU+ 
and NU+ − NT + are merged together. We use the rightmost edge of the peak in χ0 to 
determine the transition temperature. Also note the fluctuating behavior of χ2 and χp 
in the NT + phase and at the onset of N∗ (c, d), explained in section 4.3.3. <BR><BR>
T 
a) Uniaxial and biaxial order parameters q0 and q2, as defined in (3.21). <BR><BR>b) Tetrahedral order parameter Δ(3) calculated as (4.6), phase chirality calculated as <BR><BR>
22 
average parity p from (4.14) and susceptibility of tetrahedratic order parameter χ (3)
T
2 
calculated as (4.12). <BR><BR>c) Susceptibilities of uniaxial and biaxial order parameters χ0 and χ2, calculated as 

(3.44) and (3.45). <BR><BR>d) Fluctuation of energy c calculated as (3.42) and susceptibility of parity order param-<BR><BR>eter χp calculated as (4.11). <BR><BR>
χ0 
q 0 , q 2 
χp χ T (3)
2 
ab 

1 

1 

0.0007 
0.8 
0.0006 
0.8 
0.6 
0.0005 
0.4 0.6 
)(3
χ2 
Δ 
22 , p 
c 
0.0004 
0.2 0.00030.4
0 
0.0002 
-0.2 
0.2 
0.0001
-0.4 
0 
0
-0.6 
cd 

7e-05 

0.0009 30 

0.0018 
0.0008 
0.0016
6e-05 
25 
0.0007 
0.0014 
5e-05 
0.0006 
0.0012
20 
4e-05 
0.0005 
0.001 
15 
0.0004 
0.0008
3e-05 
0.0003 
0.0006
102e-05 0.0002 
0.0004 1e-05 
5
0.0001 
0.0002 
0 0 0
0 
1
Figure 4.6: Monte Carlo results for τ =1 and λ = √ . Obtained transition tempera
6 
tures: tI−NB ≈ 1.2,tNB −N∗ ≈ 1.12. Note that the NB phase exists in a very narrow 
T 
temperature range, therefore the peak in c for the I − NB transition is not discernible from the one at NB − N∗ transition.
T 
a) Uniaxial and biaxial order parameters q0 and q2, as defined in (3.21). b) Tetrahedral order parameter Δ(3) calculated as (4.6), phase chirality calculated as 
22 
average parity p from (4.14) and susceptibility of tetrahedratic order parameter χ (3)
T
2 
calculated as (4.12). <BR><BR>c) Susceptibilities of uniaxial and biaxial order parameters χ0 and χ2, calculated as 

(3.44) and (3.45). <BR><BR>d) Fluctuation of energy c calculated as (3.42) and susceptibility of parity order param-<BR><BR>eter χp calculated as (4.11). <BR><BR>
ab 

0.4 
0.2 

0.8 
1 

0.002 
0.0015 
0 
0.6 
)(3
χ 2 
22 , p
-0.2 
0.001 
-0.4 
Δ

0.4 
-0.6 
0.0005 
0.2
-0.8 -1 0 0 
cd 

8e-05 

0.0006 18 

0.012 
16
7e-05 
0.0005 
0.01 
14
6e-05 
0.0004 
0.008
12
c

5e-05 
10 
4e-05 
0.0003 
0.006 
8 
3e-05 
0.0002 
0.004
6 
2e-05 
4 
0.0001 
0.002 
1e-05 
2 
0 

0 0
0 
Figure 4.7: Monte Carlo results for τ =1 and λ =0.5. Obtained transition temperatures: tI−NU− ≈ 1.4,tNU−−NT − ≈ 1.24,tNT −−N∗ ≈ 0.94. Note the fluctuating behavior 

T 
of χ2 and χp in the NT + phase and at the onset of N∗ (c, d), explained in section 4.3.3. 
T 
a) Uniaxial and biaxial order parameters q0 and q2, as defined in (3.21). <BR><BR>b) Tetrahedral order parameter Δ(3) calculated as (4.6), phase chirality calculated as <BR><BR>
22 
average parity p from (4.14) and susceptibility of tetrahedratic order parameter χ (3)
T
2 
calculated as (4.12). <BR><BR>c) Susceptibilities of uniaxial and biaxial order parameters χ0 and χ2, calculated as 

(3.44) and (3.45). <BR><BR>d) Fluctuation of energy c calculated as (3.42) and susceptibility of parity order param-<BR><BR>eter χp calculated as (4.11). <BR><BR>
χ 0 
q 0 , q 2 
χp χ T (3)
2 
4.5 Results for τ = 28 
15 
This result corresponds to the phase diagram in Fig. 4.4b. Temperature scans, including behavior of order parameters q0, q2, p and Δ(3), energy fluctuation c and susceptibilities 
22 
χ0, χ2, χp and χ (3) are presented in figures 4.8-4.10, where transition temperatures are 
T
2 
indicated. The obtained transition temperatures have been indicated in Fig. 4.4b as 
1
points for λ =0.3, λ = √ and λ =0.5. Below will follow a general description of the 
6 
results, while detailed comments on the behavior of the measured quantities are given 
1
in the captions of figures: Fig. 4.8 for λ =0.3, Fig. 4.9 for λ = √ and Fig. 4.10 for 
6 
λ =0.5. 
1
For λ = √ the multicritical Landau point exists (Fig. 4.9), where there is a direct, 
6 
sharply first-order transition from isotropic I to the chiral phase N∗ at tI−N∗ ≈ 1.57.
T 
T 
All measured order parameters experience a discontinuous increase at this point. 
For the generic prolate (λ =0.3, Fig. 4.8) and generic oblate (λ =0.5 4.10) cases, as the temperatureislowered, asharplyfirst-orderphasetransitiontothenematictetrahedratic 
phase occurs, evidenced by a discontinuous increase in Δ(3) and q0, coupled with 
NT 22 peaks in c, χ (3) and χ0. Because of the slow dynamics of parity in lower temperatures 
T
2 
and because of the fact that biaxial ordering in the presence of tetrahedratic order is linked to parity, p and q2 experience fluctuations due to long temporal correlations in the NT phase, evidenced by the behavior of χp and χ2. The correlations continue to the N∗ phase, but wane in lower temperatures, which is visible in χp. The transition 
T 
from NT to N∗ is second-order and is evidenced by an increase in p and q2. In the case 
T 
of λ =0.5 it was possible to unequivocally identify a peak in χp, however this was not possible for λ =0.3 and an approximation is made by observing the temperature at which p begins to rise. Note that the T and NU− phases, predicted by mean-field theory 
11
(see Fig. 4.4b) to precede the NT phase for λ< √ and λ> √ respectively, exist in a 
66 
very narrow temperature range and are not observed in simulation. 
χ 0 
q 0 , q 2 
χ p χ T (3)
2 
ab 

1 

1 

0.0003 
0.8 
0.00025
0.8 
0.6 
0.0002 
0.60.4 0.00015
)(3
Δ 
22 , p 
0.2 0.4 0.0001
0 
0.2
-0.2 
5e-05 
-0.4 
0 
0 
cd 

9e-05 

0.00012 35 

0.001 
8e-05 
30
0.0001 
0.0008
7e-05 
25 
6e-05 
8e-05 
c 

0.0006
20
5e-05 
6e-05 
4e-05 
15 
0.0004 
χ2
3e-05 
4e-05 
10 
2e-05 0.0002
2e-05 1e-05 
5 
0 0 0
0 
28
Figure 4.8: Monte Carlo results for τ = and λ =0.3. Obtained transition tempera
15 
tures: tI−NT + ≈ 1.56,tNT +−N∗ ≈ 1.1. Note that the T phase, predicted by mean-field 
T 
(Fig. 4.4b) remains undetected due to very narrow temperature range. Also note the <BR><BR>fluctuating behavior of χ2 and χp in the NT + phase and at the onset of N∗ (c, d),

T 
explained in section 4.3.3. In this case tNT +−N∗ was determined from the behavior of p

T 
(b). 
a) Uniaxial and biaxial order parameters q0 and q2, as defined in (3.21). <BR><BR>
b) Tetrahedral order parameter Δ(3) calculated as (4.6), phase chirality calculated as 
22 
average parity p from (4.14) and susceptibility of tetrahedratic order parameter χ (3)
T
2 
calculated as (4.12). <BR><BR>c) Susceptibilities of uniaxial and biaxial order parameters χ0 and χ2, calculated as 

(3.44) and (3.45). <BR><BR>d) Fluctuation of energy c calculated as (3.42) and susceptibility of parity order param-<BR><BR>eter χp calculated as (4.11). <BR><BR>
χ0 q 0 , q 2 
χp χ T (3)
2 
ab 

1 

1 

0.00016 
0.00014
0.8 
0.8 
0.00012
0.6 
0.00010.4 0.6 
)(3
χ2 
Δ 
22 , p 
c 
8e-050.2 0.4 6e-05
0 
4e-05
-0.2 
0.2 
2e-05
-0.4 
0 
0
-0.6 
cd 
3.5e-05 

0.0006 18 

0.0025 
16
3e-05 
0.0005 
0.002
14 
2.5e-05 
0.0004 
12 
0.0015
2e-05 
10 
0.0003 
8
1.5e-05 
0.001 
0.0002 6
1e-05 
4 
0.0005
0.0001
5e-06 
2 
0 

0 0
0 
28 1
Figure 4.9: Monte Carlo results for τ = and λ = √ . The multicritical point is found 
15 
6 
for tI−N∗ ≈ 1.57. 
T 
a) Uniaxial and biaxial order parameters q0 and q2, as defined in (3.21). b) Tetrahedral order parameter Δ(3) calculated as (4.6), phase chirality calculated as 
22 
average parity p from (4.14) and susceptibility of tetrahedratic order parameter χ (3)
T
2 
calculated as (4.12). <BR><BR>c) Susceptibilities of uniaxial and biaxial order parameters χ0 and χ2, calculated as 

(3.44) and (3.45). <BR><BR>d) Fluctuation of energy c calculated as (3.42) and susceptibility of parity order param-<BR><BR>eter χp calculated as (4.11). <BR><BR>
ab 


1 

0.00014 
0.4 
0.00012 
0.2 
0.8 
0.0001
0 0.6 
8e-05 
χ 0 q 0 , q 2 
χ p χ T (3)
2 
-0.2 
)(3
Δ 
22 , p 
-0.4 
6e-050.4 
-0.6 
4e-05 
0.2
-0.8 
2e-05 
-1 
0 
0 
cd 

3.5e-05 

0.0002 16 

0.0025 
14
3e-05 
0.002 
0.00015 12
2.5e-05 
c

10 
0.0015
2e-05 
0.0001 
8 
1.5e-05 
0.0016 
χ2 
1e-05 5e-05 
4 
0.0005 
5e-06 2 
0 
0 
0 
28
Figure 4.10: Monte Carlo results for τ = and λ =0.5. Obtained transition temper
15 
atures: tI−NT − ≈ 1.7,tNT −−N∗ ≈ 1.21. Note that the T phase, predicted by mean-field 
T 
(Fig. 4.4b) remains undetected due to very narrow temperature range. Also note the fluctuating behavior of χ2 and χp in the NT − phase and at the onset of N∗ (c, d),
T 
explained in section 4.3.3. <BR><BR>a) Uniaxial and biaxial order parameters q0 and q2, as defined in (3.21). <BR><BR>
b) Tetrahedral order parameter Δ(3) calculated as (4.6), phase chirality calculated as 
22 
average parity p from (4.14) and susceptibility of tetrahedratic order parameter χ (3)
T
2 
calculated as (4.12). <BR><BR>c) Susceptibilities of uniaxial and biaxial order parameters χ0 and χ2, calculated as 

(3.44) and (3.45). <BR><BR>d) Fluctuation of energy c calculated as (3.42) and susceptibility of parity order param-<BR><BR>eter χp calculated as (4.11). <BR><BR>
0 

4.6 Discussion and summary 
Two classes of phase diagrams for the model (4.1) have been studied by Monte Carlo simulation, providing the dependence on temperature and model parameters λ and τ of order parameters, fluctuation of energy and susceptibilities (measures of fluctuation) of order parameters for typical prolate and oblate values (λ =0.3 and λ =0.5), as well as the highly-symmetric Landau point (self-dual point of the quadrupolar biaxial model 
1
(3.1)), λ = √ . The two studied values of τ pertain to the nematic-dominated phase 
6 28
diagram topology for τ � 1.54 (τ =1) and the highly-symmetric value τ = , for which 
15 
the multicritical Landau point is observed. 
Apart from the phases existing in l =2 model of biaxial nematics (3.1), two phases with tetrahedratic ordering have been identified: nematic tetrahedratic NT and chiral nematic N∗ . The T phase was not accessible by simulation in the considered cases. For τ =1
T 
28
it is not present on the phase diagram, while for τ = its temperature range is too 
15 
narrow to allow identification of the phase. 
In the NT phase spontaneous breaking of tetrahedratic symmetry was observed by measuring the tetrahedratic order parameter Δ(3). The phase transition to NT also is seen to 
22 
produce pronounced peaks in the fluctuation of energy c and susceptibility of the tetrahedratic order parameter χ (3) . The N∗ phase was identified by observing the spontaneous 
T
T
2 
breaking of chirality, measured in terms of average parity p. The chiral phase is observed in all cases, however the transition temperature is highest for a highly biaxial molecular 
1
quadrupolar moment (λ ≈ √ ) and high tetrahedratic coupling (through τ).
6 
The results need not apply exclusively to the speculated complexes of embrace-locked bent molecules, as in Fig. 4.1, and can be linked to cases where the leading terms in the dispersion interaction expansion are the quadrupolar and octupolar couplings. To provide a general correspondence, in the future also tetrahedratic coupling terms where the octupolar tensor’s twofold axes are rotated with respect to the twofold axes of the quadrupole need to be considered. However, the observed homogeneous chirality without the observation of twisted domains or macroscopic states seems unlikely in real compounds. (4.1) does not account for twisted structures because of the lack of spatially<BR>dependent coupling (e.g. coupling of the quadrupolar or octupolar moments to lattice vectors). This type of coupling is accounted for when, in addition to E> 0, τ> 0 and λ> 0 in the general model (2.11) κ> 0, the case studied in Chapter 5. 
The results also indicate that the investigation of spontaneous chirality breaking in liquid crystalscallsfornewsimulationmethods. Inourexperience,brute-forceincreasingofthe cycle count in a standard Metropolis algorithm provides less-than-satisfying output when compared to the increasing time and computer resource consumption when a discrete degree of freedom, such as parity, is taken into account. Therefore in the future methods such as parallel tempering should be considered [123]. 
Chapter 5 


Spatially modulated structures 
To conclude this thesis with a promising, open problem, we explore a proposed extension to the model (4.1) which, as we preview in this short chapter, produces spatially inhomo<BR>geneous structures and ambidextrous chirality. Such an effect should be accounted for, considering the findings involving bent-core molecules [38–40, 43, 55] and flexible dimers [45–50]. 
5.1 Intermolecular vector coupling 
The inclusion of uniaxial (D∞h), biaxial (D2h) and tetrahedral (Td) symmetries in the model (4.1) permits one to construct an antisymmetric tensor using the symmetrized irreducible tensors: 
√ 
(2) (2) (3)
ˆ
bαn · (ˆbβn × ˆbγn)=2 2 T0,xµ T2,yν T2,µνz, (5.1) (x,y,z)∈π{α,β,γ} 
ˆˆˆ
where ˆbxn = aˆn,byn = bn,bzn = cˆn are the molecular basis vectors for site n and the summation runs over cyclic permutations of {α, β, γ}, with Einstein summation convention assumed for the remaining indices. The lowest non-trivial term including this tensor 
involves coupling with the intermolecular (or inter-site) versor rˆij = rij (nrij = nrj − nri),
rij 
rˆij = −rˆji and one of the molecular moments present in the interaction (4.1). Consider<BR>ing the quadrupolar moment, one obtains [119]: 
Vc(Ωˆi, Ωˆj)= κE [ˆbαi · (ˆbβi × ˆbγi)Qαν (Ωˆi)Qβν(Ωˆj) 
(5.2) 
− ˆbαj · (ˆbβj × ˆbγj)Qαν (Ωˆj)Qβν (Ωˆi) (ˆrij)γ, 
101 

where α, β, γ, ν = {x, y, z}, Qρδ(ˆΩ) 
Ω) = Q(ˆand Einstein summation convention is 
ρδ 
assumed. A similar term exists for the octupolar tensor: 
(3) (3)
VcT (Ωˆi, Ωˆj)= κT E [ˆbαi · (ˆbβi × ˆbγi)T2 ανµ(Ωˆi)T2 βνµ(Ωˆj) 
(5.3)
(3) (3)
− ˆbαj · (ˆbβj × ˆbγj)T2 ανµ(Ωˆj)T2 βνµ(Ωˆi) (ˆrij)γ, 
(3) (3)
(ˆ(ˆ
where T Ω) = TΩ) . Terms (5.2) and (5.3) are invariant with respect to 
2 αρδ2 ρδ
simultaneous inversion of both parities and an exchange of sites, i.e.: 
V (T )(pi, Ωˆi' ,pjΩˆj' )= V (T )(−pj, Ωˆj' , −pi, Ωˆi' ). (5.4)
cc 

5.2 Consequences of intermolecular vector coupling 
We study a model defined by the pair interaction potential constructed by supplementing the potential (4.1) with the quadrupolar variety of the lattice-coupling potential (5.2): 
ˆˆˆ
Vt(Ωˆi, Ωj)= VT (Ωˆi, Ωj)+ Vc(Ωˆi, Ωj). (5.5) 
We assume κ ≥ 0. The Hamiltonian is constructed on a lattice with nearest-neighbors coupling: 
N
1 Vt(ˆˆ
H =Ωi, Ωj). (5.6)
2 
(i,j) 
This model corresponds to the general model (2.11) for the case when the parameters λ, E, τ and κ are non-zero. We do not specify which type of lattice is assumed, as three types are used: an open chain, a two-dimensional square lattice and a three-dimensional simple cubic lattice. 
We will now consider the minimum of (5.5) for two isolated neighboring particles to show that it favors a twisted state for non-zero κ. The simplest parametrization of (5.5) 
ˆΩˆ' Ω '
is in relative Euler angles, i.e. Vt(Ωˆi, Ωj) ≡ Vt(pi,pj, ), where ˆ= {αij,βij,γij}
ijij 
defines the (proper) rotation, which carries Ωˆ' j into Ωˆ' i. Since i and j are fixed, for simplicity we assume ˆΩ= {α, β, γ}. For the sake of the argument we place the 
Ω ' ≡ ˆ
ij 
particle i at the origin and the intermolecular vector, which points at the particle j, is parametrized in polar coordinates (θ, φ): rˆij = {sin(θ) cos(φ), sin(θ) sin(φ), cos(θ)}. In this parametrization, (5.5) reads: 
Vt(pi,pj, Ωˆ' )= −E[fQQ(Ωˆ' )+ τpipjfTT (Ωˆ' )+ κ(pi + pj)fr(Ωˆ' , θ, φ)], (5.7) 
where: 
√
(2) (2) (2) (2)
Ω ' )=Δ(Ωˆ' )+ 2λ[Δ(Ωˆ' )+Δ(Ωˆ' )] + 2λ2Δ(Ωˆ' ), (5.8)
fQQ(ˆ00 2002 22 
(3)
fTT (Ωˆ' )=Δ22 (Ωˆ' ), (5.9) 
√
1 
fr(Ωˆ' , θ, φ)= {λ sin(2α cos(θ)[λ cos(2γ)(3 + cos(2β))+ 6sin2(β)]+
2
√ 
+ λ sin(2γ)[sin(β) sin(θ)(2λ cos(α + φ) + 6 cos(α + φ))+ 
(5.10)
1 
+4λ cos(2α) cos(β) cos(θ)] + sin(2β) sin(θ)[sin(α − φ)× 
2 
√√ × (6λ cos(2γ) − 3) − λ sin(α + φ)( 6 − 2λ cos(2γ))]}, 
(3) (2) (2)
Δ(Ωˆ' ), Δ(Ωˆ' ) and Δ(Ωˆ' ) are the symmetry-adapted Wigner functions (2.6), dis<BR>
2200 22 
cussed in section 2.1. 
For κ =0 (and E> 0,τ > 0), the minimum of (5.7) is: 
ˆ
minVt(pi,pj, Ω ' )= −E(1 + 2λ2 + pipjτ), (5.11) 
{Ωˆ' } 
which is satisfied when pi = pj and the tetrahedral and quadrupolar twofold axes are aligned in parallel. For κ �
=0, a competing minimum becomes global, which is achieved when the minor molecular axes aˆi and aˆj align in parallel along the vector rˆij, while the perpendicular molecular axes {ˆbj,cˆj} are rotated along aˆj by an angle β0ij, relative to {ˆbi,cˆi}. The rotation is clockwise if pi = pj =1 and counter-clockwise if pi = pj = −1. The relative angle is given by: 
κΛ(pi + pj)
tan(2β0ij)= , (5.12)
Λ+4τ 
where: 
√ 
Λ=3+2λ( 6+ λ). (5.13) 
(5.12) can give a general idea on the macroscopic implications of the ground state, when expanded in small κ. Upon ensemble-averaging the leading term, we obtain: 
β0ij ≈ κΛ p. (5.14)
Λ+4τ 
It is now evident that the N∗ phase, i.e. when p �0, is unstable with respect to 
= 
T 
spontaneous formation of twisted states and that the predictions of chapter 4 are valid only in the κ =0 limit. In the above approximation the twisted states have pitch of the 
π
orderof .Inotherphases,where p =0, from(5.12)wecanconcludethatlocaltwisted 
β0ij 
domains of opposite handedness will form in equal abundance, leading to ambidextrous chirality. 

�2012 by the American Physical Society.) 
5.2.1 One-dimensional chain 
The prediction (5.14) is very well supported by Monte Carlo simulations on a one<BR>dimensional open chain of molecules. We consider 32 sites aligned along the direction xˆwith free boundary conditions in all directions. The degrees of freedom are in no way limited with respect to the three-dimensional case, thus the molecules are free to rotate, as in previously considered cases. 
We perform Monte Carlo simulations with the pair interaction potential (5.5) for τ = 1,λ =0.3 and κ =1. Considering homogeneous parity, p =1, the estimation (5.14) gives β0ij ≈ 31◦ for mean relative clockwise twist about the molecular aˆaxes and a mean 
π
pitch ≈ 5.8. As can be viewed in Fig. 5.2a, an equilibrium state after a 60000-cycle 
β0ij 
simulation for temperature t =0.2, starting from an ordered state with homogeneous parity p =1, confirms these predictions. The pitch is roughly equal to 6, while the mean relative twist angle, averaged over the chain in the pictured configuration, is equal to β0ij ≈ 30.8◦. On the other hand, simulation for a higher temperature shows domains of opposite chirality with small but discernible twists, which confirms our predictions. (See Fig. 5.2b.) However, it should be stressed that although overall parity p decreases as temperature is raised, there is no phase transition observed in the one-dimensional case, as it is known that for discrete symmetry groups no phase transitions exist for a chain with nearest-neighbors interaction in finite temperature. 

Figure 5.2: Result of a 60000-cycle Monte Carlo simulation of (5.5) on an open-ended one-dimensional chain of length 32 for t =0.2 (a) and t =0.5 (b) with τ =1,λ =0.3 and κ =1. The pictured rectangular cuboids represent the molecular quadrupolar moments Q(ˆ
Ω), while the color stands for parity: p = +1 (red) and p = −1 (green). a) Global twisted state in the homochiral phase. The pitch of the global twist is estimated as ∼ 6, while the mean relative angle is β0ij ≈ 30.8◦ (averaged over the pictured configuration), which is in excellent agreement with the prediction (5.14). b) Ambidextrous chiral domains in a disordered phase. 


5.3 Two-dimensional lattice 
Let us constrain the model (5.5) to a two-dimensional square lattice in the (ˆx, yˆ) plane and let the boundary conditions in xˆand yˆdirections be periodic, while keeping the boundary conditions free along zˆ. The molecular degrees of freedom remain as in the three-dimensional case, thus we are left with a free-standing single-layered film. For such a system the intermolecular vectors are constrained to four directions, i.e. rˆij = {1, 0}, {0, 1}, {−1, 0}, {0, −1}. 
Consider now that two molecules, 1 and 2, of identical parity p1 = p2 = +1, occupy neighboring sites at positions (x1,y1) = (0, 0) and (x2,y2) = (1, 0) (Fig. 5.3a). As discussed in the previous section, the ground state of such configuration is achieved when both molecules 1 and 2 align with aˆ1 aˆ2 rˆ12 = {1, 0} and molecule 2 is tilted clockwise with respect to 2 around aˆ2 (or vice-versa). Now consider a third molecule, 3, p3 = p1 = p2 = +1, located at (x3,y3) = (0, 1). If 1 and 3 were treated in isolation (Fig. 5.3b), the ground state would be achieved analogously, by aligning aˆ3 aˆ1 rˆ13 and tilting 3 clockwise with respect to 1 along aˆ3 (or vice-versa). Now notice that the two pairwise ground states for 1 and 2 and for 1 and 3 cannot be achieved simultaneously (Fig. 5.3c). Adding more molecules multiplies the number of conflicting conditions. Thus, the system is frustrated. The question of how the system relaxes this frustration is difficult to answer on analytical grounds and in the present we resort to simulation in search for an answer. 
5.3.1 Simulations and results 
We perform Monte Carlo simulations of the model (5.5) on a square 32 × 32 lattice with periodic boundary conditions in the lattice plane and free boundary conditions in the perpendicular direction. We take λ =0.3 and τ = κ =1. Starting from an initial, 

Q(ˆ
Ω). a) Ground state of two isolated molecules placed on neighboring lattice sites along xˆ. b) Ground state of two isolated molecules placed on neighboring lattice sites along yˆ. c) The two pairwise ground states a and b conflict each other. 
completely ordered and righthanded (p = +1) configuration, we run a simulation of 60000 cycles for the lowest temperature (t0 =0.2) and pass the final state as the initial state to the next temperature and iterate until t =1.05 is reached with an increment of δt =0.05. To compare the scale of the number of Monte Carlo cycles taken into account, consider that it amounts to average ∼ 60 lattice sweeps per particle on the 32×32 lattice, which is of similar order as ∼ 50 − 100 lattice sweeps per particle on the 16 × 16 × 16 lattices for 200000 − 500000 cycles considered in the previous chapter. 
In two dimensions we do not expect any phase transition associated with the continuous (rotational) degrees of freedom. However, as discussed in chapter 4, the molecular parity, being a discrete, two-state degree of freedom, is an Ising-like variable. Therefore it is reasonable to expect a phase transition with spontaneous parity (chirality) breaking. Certainly, chiral ordering will consequently induce some level of biaxial and tetrahedral ordering, but this ordering is not due to spontaneous breaking of SO(3) symmetry. We monitor this phase transition by the fluctuation of energy c (3.42) and the susceptibility of the parity order parameter χp (4.11). Average parity is computed by direct ensemble averaging of the molecular parity p. All other quantities are computed as they have been in previous chapters. 
The results indicate that, indeed, a phase transition between the disordered and chiral states exists in two dimensions, at temperature tc ≈ 0.9 ∼ 0.95, characterized by a sharp drop in phase chirality p and peaks in c and χp as tc is approached from below (Fig. 5.4). This is an expected result, however what is more curious is the apparent existence of a phase transition within the chiral phase. It is visible as a minor peak in c (Fig. 5.4b) and a pronounced peak in the susceptibility of the uniaxial order parameter χ0 (3.44), as well as an increase in the uniaxial parameter q0 (3.21) as the temperature is raised 
p, q 0
χp 
ab 
5e-05 12 
1 
0.005 
10
4e-05 
0.8 
0.004 

3e-05 

0.002 
4 
8 
χ 0 
c 
0.003 
6
2e-05 
1e-050.2 
0.001 
2 
0 0 0 
Figure 5.4: Several Monte Carlo quantities obtained from a simulation of model (5.5) on <BR><BR>a 32 × 32 lattice for λ =0.3 and τ = κ =1. 
a) Phase chirality measured by average parity p, the primary uniaxial order parameter 
q0 (3.21) and susceptibility of the uniaxial order parameter χ0 (3.44). <BR><BR>b) Fluctuation of energy c (3.42) and susceptibility of the chirality order parameter χp 
(4.11). <BR><BR>
(Fig. 5.4a), which indicates that some kind of reordering is taking place, approximately at the temperature tr ≈ 0.55. 
Investigation of the final lattice states reveals that for t<tr bent-wavefront-cholesteric structures exist, where the cholesteric-like wavefronts bend at right angles and do not form any defects, i.e. do not split or terminate (Fig. 5.5d). This structure apparently relaxes the frustration of the ground state, discussed in the previous section. We have found these states to be remarkably stable, as they form rapidly from the initial ordered state and throughout the simulation no qualitative changes are observed. For tc >t>tr two-dimensional cholesteric structures are found, where the molecules are aligned with their molecular aˆaxes in parallel to the diagonal of the lattice and rotate along aˆ(Fig. 5.5b). The pitch of such structure appears to agree with the prediction (5.14) in the xˆ
1π
and yˆdirections, while along the diagonal it can be estimated as √ , which for the 
2 (β0ij )
given parameters predicts a pitch of 4. Investigation of e.g. Fig. 5.5b shows that this is indeed the case. Above tc, the disordered state is found with many small domains of opposite chirality (Fig. 5.5a). The ordered structures bear uncanny resemblance to the freeze-fracture imaging results obtained recently for the twist-bend nematic phase [50]. 

Ω), while the color represents parity: p = +1 (red) and p = −1 (green). a) t =0.95: Ambidextrous chiral domains in the disordered state existing for t>tc. b)t =0.8: Atwo-dimensionalcholestericstatewithawavevectorparalleltothediagonal of the lattice. The molecules are aligned with their molecular aˆaxes in parallel, along which they rotate. c) t =0.55: Reordering transition in the chiral phase. At t = tr ≈ 0.55 cholesteric states found for tc >t>tr coexist with the bent-wavefront-cholesteric structures encountered for t<tr. d) t =0.25: Stable bent-wavefront-cholesteric structures existing for t<tr. 


5.4 Three-dimensional lattice 
The results for the two-dimensional version of model (5.5) raise the question whether similar bent-wavefront structures, or perhaps even more complicated ones, will be observed in three dimensions. Clearly, the frustration of the ground state argued for the two-dimensional case should be more manifest in three dimensions. Our results in this respect are merely a teaser of what is possible in these systems, therefore we keep this section short and skip to presentation of the results. 
The simulations were carried out for a three-dimensional 16×16×16 simple cubic lattice with periodic boundary conditions in all directions. We proceed in a similar manner as in the two-dimensional case, i.e. by first providing an initial, completely ordered state with p = +1 and passing the final state of each simulation for a given temperature as the initial state of the subsequent simulation for a higher temperature. The considered parameters are λ =0.3 and κ = τ =1, while the considered cycle count is 80000. 
These preliminary results indicate that also in the case of a three-dimensional lattice, complex, frustrated supramolecular states may exist. They have been found in the lower temperature regime, e.g. t =0.2 and t =0.4 (Fig. 5.6a and b), however for higher temperatures a periodic phase is also encountered, e.g. for t =1.0 and t =1.2 (Fig. 5.6c and d). As we see in Fig. 5.6c in particular, the molecules twist in a helix approximately along the diagonal and the pitch can be assessed by the modulating color to span only a few molecules. This kind of ordering is typical of the twist-bend nematic phase [42], while the surfaces of the simulation box in Fig. 5.6c and d can be, as in the 2D case, compared to the freeze-fracture imaging results obtained this year by the group at the University of Colorado in Boulder [50]. According to this account, the striped textures of the frozen samples do not indicate layering, but reflect the orientation of bent molecules forming helices. This is also visible in Fig. 5.6c and d. 

Ω), while the color represents parity: p = +1 (red) and p = −1 (green) with hue modulated by orientation of the molecular cˆaxes. a,b) t =0.2 and t =0.4: Stable complex supramolecular structures. c,d) t =1.0 and t =1.2: Stable layered phases. 

5.5 Discussion 
The preliminary results presented in this chapter indicate that the complex effects due to spontaneous breaking of chirality, encountered e.g. in bent-core and flexible dimer systems can be accounted for in a microscopic dispersion model with a coupling to the intermolecular vectors, in the form such as (5.2). So far, such a model, in which twisted states and supramolecular structures form due to spontaneous breaking of chirality, has been missing. Subsequent detailed analysis needs to be undertaken to study the phase transitions and remaining unidentified structures, which can exist in other phases, in particular in the three-dimensional case. Finally, because geometrical frustration and chirality are known to lead to the emergence of blue phases, it is important to establish whether a link between the model (5.5) to those structures exists. 
Chapter 6 


Summary and conclusions 

In the present thesis we have analyzed some of the cases described by the general model 
(2.11) and have shown that it does reproduce many of the remarkable effects observed in bent-shaped liquid crystal molecules (bent-core, flexible dimer and ferrocene mesogens), such as field-stabilized biaxiality, homogeneous chirality, tetrahedratic order and twisted states of homogeneous and ambidextrous chirality. The results have been summarized in detail at the conclusions of the respective chapters, 3, 4, and 5. Herein we present a general comment on the results when viewed as a whole. 
In general, the model studied in the present thesis, (2.11), produces many of the above phenomena, such as field-stabilized biaxiality (as seen Chapter 3, in section 3.6.2, in particular Fig. 3.16), homogeneous chirality and tetrahedratic phases (Chapter 4, in sections 4.4 and 4.5), twisted helical states (in particular the three-dimensional result in section 5.4 in Chapter 5) and ambidextrous chirality (corroborated by 1D and 2D results in sections 5.2.1 and 5.3 in Chapter 5), as well as effects present in classic nematogens, such as the critical and tricritical points of the isotropic-uniaxial nematic transition (studied in detail in section 3.6.1 in Chapter 3). The stable biaxial nematic phase, which is present for no external field (see section 3.5) and also for τ> 0 (in particular see Fig. 4.6), has provedelusive in experimentand is widely believed to be unstableand is subject to constant experimental research at the time of completion of this thesis. Some cases of the model (2.11) however produce novel results, such as the field-induced decrease in 
1
the isotropic-biaxial nematic transition temperature for the model (3.14) at λ = √ (see the mean-field phase diagram in Fig. 3.18 and order parameter comparison in Figs. 3.19 and 3.20) and the apparent emergence of complex, frustrated supramolecular structures produced by the model (5.5), presented in Fig. 5.6 in Chapter 5, which bear resemblance to the purported twist-bend nematic phase and blue phases. 
111 

Out of the vast number of possible cases we have investigated merely a few. Before moving on to other ones than those considered in the present thesis, further investigation into the structures produced by the model (5.5), which we begin to explore in Chapter 5, is desirable. Other cases of great interest have been omitted here, which includes prominently the case of E> 0,λ > 0,τ > 0 and h> 0 in (2.11) (other parameters being equal to zero), in which the external field is expected to induce chirality in the nematic tetrahedratic (NT ) phase. Furthermore, if one additionally considers κ> 0, if field-induced chirality is indeed found, a possibility of field-triggered modulated structures opens, which clearly deserves attention for practical purposes. 
Thanks 
I would like to thank my supervisor, Prof. Lech Longa, for his teachings, support, relentlessness, uplifting optimism and guidance throughout my studies. 
I thank Dr. Tomasz Wydro for his indispensable help, company and many inspiring discussions. 


Bibliography 
[1] Longa, Pająk, and Wydro. Chiral symmetry breaking in bent-core liquid crystals. Physical Review E, 79(4):040701, 2009. 
[2] de Gennes and Prost. The physics of liquid crystals. Clarendon Press (Oxford and New York), 1993. 
[3] Neto and Salinas. 	The Physics of Lyotropic Liquid Crystals: Phase Transitions and Structural Properties, volume 62. OUP Oxford, 2005. 
[4] TolédanoandNeto. Phase transitions in complex fluids. WorldScientificPublishing Company, 1998. 
[5] Reinitzer. Beiträge zur kenntniss des cholesterins. Monatshefte für Chemie/Chemical Monthly, 9(1):421–441, 1888. 
[6] Lehmann. Über fliessende krystalle. Z. phys. Chem, 4:462–472, 1889. 
[7] Lehmann. Structure, system and magnetic behaviour of liquid crystals and their miscibility in solid crystals. Annalen der Physik (series 4), 2:649–705, 1900. 
[8] Priestley and Wojtowicz. Introduction to liquid crystals. Plenum Publishing Corporation, 1975. 
[9] Chandrasekhar. Liquid crystals. Cambridge University Press, 1992. 
[10] Vertogen and de Jeu. 	Thermotropic liquid crystals, fundamentals, volume 45. Springer-Verlag Berlin, 1988. 
[11] Heilmeier. Liquid-crystal display devices. Scientific American, 222:100–106, 1970. 
[12] Bahadur. Liquid crystal displays. Molecular Crystals and Liquid Crystals, 109(1): 3–93, 1984. 
[13] Chigrinov. Liquid crystal devices: physics and applications. Artech House Publishers, 1999. 
115 

[14] Boccara. 	Violation of rotational invariance and mesomorphic phase transitions characterized by an order parameter. Annals of Physics, 76(1):72–79, 1973. 
[15] Tinh, Destrade, and Gasparoux. Nematic disc-like liquid crystals. Physics Letters A, 72(3):251–254, 1979. 
[16] de Vries. Rotatory power and other optical properties of certain liquid crystals. Acta Crystallographica, 4(3):219–226, 1951. 
[17] Meiboom, Sethna, Anderson, and Brinkman. 	Theory of the blue phase of cholesteric liquid crystals. Physical Review Letters, 46(18):1216–1219, 1981. 
[18] Kikuchi, Yokota, Hisakado, Yang, and Kajiyama. Polymer-stabilized liquid crystal blue phases. Nature Materials, 1(1):64–68, 2002. 
[19] Yu and Saupe. 	Observation of a biaxial nematic phase in potassium laurate-1decanol-water mixtures. Physical Review Letters, 45(12):1000–1003, 1980. 
[20] Severing and Saalwächter. 	Biaxial nematic phase in a thermotropic liquidcrystalline side-chain polymer. Physical review letters, 92(12):125501, 2004. 
[21] Freiser. Ordered states of a nematic liquid. Physical Review Letters, 24(19):1041– 1043, 1970. 
[22] Hessel, Herr, and Finkelmann. Synthesis and characterization of biaxial nematic sidechainpolymerswithlaterallyattachedmesogenicgroups. Die Makromolekulare Chemie, 188(7):1597–1611, 1987. 
[23] Chandrasekhar, Sadashiva, Ratna, and Raja. A biaxial nematic liquid crystal. Pramana-Journal of Physics, 30(5):L491–L494, 1988. 
[24] Praefcke, Kohne, Singer, Demus, Pelzl, and Diele. Thermotropic biaxial nematic phases with negative optical character [1]. Liquid Crystals, 7(4):589–594, 1990. 
[25] Fan, Fletcher, Gündogan, Heaton, Kothe, Luckhurst, and Praefcke. The symmetry ofthenematicphaseofathermotropicliquidcrystal: biaxialoruniaxial? Chemical physics letters, 204(5):517–523, 1993. 
[26] Li, Percec, Rosenblatt, and Lavrentovich. Biaxiality in a cyclic thermotropic nematic liquid crystal. EPL (Europhysics Letters), 25(3):199, 1994. 
[27] Madsen, Dingemans, Nakata, and Samulski. Thermotropic biaxial nematic liquid crystals. Physical review letters, 92(14):145505, 2004. 
[28] Merkel, Kocot, Vij, Korlacki, Mehl, and Meyer. 	Thermotropic biaxial nematic phase in liquid crystalline organo-siloxane tetrapodes. Physical review letters, 93 (23):237801, 2004. 
[29] Acharya, Primak, and Kumar. Biaxial nematic phase in bent-core thermotropic mesogens. Physical review letters, 92(14):145506, 2004. 
[30] Prasad, Kang, Suresh, Joshi, Wang, and Kumar. Thermotropic uniaxial and biaxial nematic and smectic phases in bent-core mesogens. Journal of the American Chemical Society, 127(49):17224–17227, 2005. 
[31] Lee, Lim, Kim, and Jin. Dynamics of electro-optical switching processes in surface stabilizedbiaxialnematicphasefoundinbent-coreliquidcrystal. Journal of applied physics, 101(3):034105–034105, 2007. 
[32] Teixeira, Masters, and Mulder. Biaxial nematic order in the hard-boomerang fluid. Molecular Crystals and Liquid Crystals, 323(1):167–189, 1998. 
[33] Galerne. Comment on “thermotropic biaxial nematic liquid crystals”. 	Physical review letters, 96(21):219803, 2006. 
[34] Vanakaras and Photinos. Thermotropic biaxial nematic liquid crystals: Spontaneous or field stabilized? The Journal of chemical physics, 128:154512, 2008. 
[35] Matsuzaki and Matsunaga. 	New mesogenic compounds with unconventional molecular structures 1, 2-phenylene and 2, 3-naphthylene bis [4-(4alkoxyphenyliminomethyl) benzoates] and related compounds. Liquid Crystals, 14(1):105–120, 1993. 
[36] Takezoe and Takanishi. Bent-core liquid crystals: their mysterious and attractive world. Japanese journal of applied physics, 45:597, 2006. 
[37] Francescangeli, Vita, Ferrero, Dingemans, and Samulski. Cybotaxis dominates the nematic phase of bent-core mesogens: a small-angle diffuse x-ray diffraction study. Soft Matter, 7(3):895–901, 2011. 
[38] Link, Natale, Shao, Maclennan, Clark, Körblova, and Walba. Spontaneous formation of macroscopic chiral domains in a fluid smectic phase of achiral molecules. Science, 278(5345):1924–1927, 1997. 
[39] Xu, Selinger, Selinger, and Shashidhar. Monte carlo simulation of liquid-crystal alignment and chiral symmetry-breaking. The Journal of Chemical Physics, 115: 4333, 2001. 
[40] Yan, Hixson, and Earl. Self-assembled chiral superstructures composed of rigid achiral molecules and molecular scale chiral induction by dopants. Physical review letters, 101(15):157801, 2008. 
[41] Reddy and Tschierske. Bent-core liquid crystals: polar order, superstructural chirality and spontaneous desymmetrisation in soft matter systems. Journal of materials chemistry, 16(10):907–961, 2006. 
[42] Dozov. 	On the spontaneous symmetry breaking in the mesophases of achiral banana-shaped molecules. EPL (Europhysics Letters), 56(2):247, 2001. 
[43] Peroukidis, Vanakaras, and Photinos. Molecular simulation of hierarchical structures in bent-core nematic liquid crystals. Physical Review E, 84(1):010702, 2011. 
[44] Kim, Kadkin, Kim, and Choi. Tetrahedratic mesophases, ambidextrous chiral domains and helical superstructures produced by achiral 1, 1’-disubstituted ferrocene derivatives. European Journal of Inorganic Chemistry, 2011(19):2933–2941, 2011. 
[45] Görtz, Southern, Roberts, Gleeson, and Goodby. Unusual properties of a bent-core liquid-crystalline fluid. Soft Matter, 5(2):463–471, 2009. 
[46] Panov, Nagaraj, Vij, Panarin, Kohlmeier, Tamba, Lewis, and Mehl. Spontaneous periodic deformations in nonchiral planar-aligned bimesogens with a nematicnematic transition and a negative elastic constant. Physical review letters, 105 (16):167801, 2010. 
[47] Cestari, Diez-Berart, Dunmur, Ferrarini, de La Fuente, Jackson, Lopez, Luckhurst, Perez-Jubindo, Richardson, et al. Phase behavior and properties of the liquidcrystal dimer 1”, 7”-bis (4-cyanobiphenyl-4’-yl) heptane: A twist-bend nematic liquid crystal. Physical Review E, 84(3):031704, 2011. 
[48] Henderson and Imrie. Methylene-linked liquid crystal dimers and the twist-bend nematic phase. Liquid Crystals, 38(11-12):1407–1414, 2011. 
[49] Beguin, Emsley, Lelli, Lesage, Luckhurst, Timimi, and Zimmermann. The chirality of a twist–bend nematic phase identified by nmr spectroscopy. The Journal of Physical Chemistry B, 116(27):7940–7951, 2012. 
[50] Chen, Porada, Hooper, Klittnick, Shen, Tuchband, Korblova, Bedrov, Walba, Glaser, et al. Chiral heliconical ground state of nanoscale pitch in a nematic liquid crystal of achiral molecular dimers. Proceedings of the National Academy of Sciences, 110(40):15931–15936, 2013. 
[51] Lubensky and Radzihovsky. Theory of bent-core liquid-crystal phases and phase transitions. Physical Review E, 66(3):031704, 2002. 
[52] Mettout.Macroscopicandmolecularsymmetriesofunconventionalnematicphases. Physical Review E, 74(4):041701, 2006. 
[53] Wiant, Neupane, Sharma, Gleeson, Sprunt, Jakli, Pradhan, and Iannacchione. Observation of a possible tetrahedratic phase in a bent-core liquid crystal. Physical Review E, 77(6):061701, 2008. 
[54] Weissflog, Schröder, Diele, and Pelzl. Field-induced formation of the polar smcp phase above the smcp–isotropic transition. Advanced Materials, 15(7-8):630–633, 2003. 
[55] Brand,Pleiner,andCladis.Tetrahedraticcross-couplings:novelphysicsforbanana liquid crystals. Physica A: Statistical Mechanics and its Applications, 351(2):189– 197, 2005. 
[56] Eisenshitz and London. Zeitschrift für Physik, 60:491, 1930. 
[57] Buckingham, Fowler, and Hutson. Theoretical studies of van der waals molecules and intermolecular forces. Chemical Reviews, 88(6):963–988, 1988. 
[58] Buckingham. Theory of long-range dispersion forces. Discuss. Faraday Soc., 40: 232–238, 1965. 
[59] Maier and Saupe. Eine einfache molekulare Theorie des nematischen kristallinflüssigen Zustandes. Zeitschrift Naturforschung Teil A, 13:564, July 1958. 
[60] Maier and Saupe. 	Eine einfache molekularstatistische Theorie der nematischen kristallinflüssigen Phase. Teil I. Zeitschrift Naturforschung Teil A, 14:882, October 1959. 
[61] Maier and Saupe. 	Eineeinfache molekularstatistische Theorie der nematischen kristallinflüsign Phase. Teil II. Zeitschrift Naturforschung Teil A, 15:287, April 1960. 
[62] Onsager. The effects of shape on the interaction of colloidal particles. Annals of the New York Academy of Sciences, 51(4):627–659, 1949. 
[63] Luckhurst, Zannoni, Nordio, and Segre. 	A molecular field theory for uniaxial nematicliquidcrystalsformedbynon-cylindricallysymmetricmolecules. Molecular Physics, 30(5):1345–1358, 1975. 
[64] Luckhurst and Romano. Computer simulation studies of anisotropic systems: Ii. uniaxial and biaxial nematics formed by non-cylindrically symmetric molecules. Molecular Physics, 40(1):129–139, 1980. 
[65] Lebwohl and Lasher. 	Nematic-liquid-crystal order—a monte carlo calculation. Physical Review A, 6(1):426, 1972. 
[66] Fabbri and Zannoni. 	A monte carlo investigation of the lebwohl-lasher lattice model in the vicinity of its orientational phase transition. Molecular Physics, 58 (4):763–788, 1986. 
[67] Straley. 	Ordered phases of a liquid of biaxial particles. Physical Review A, 10: 1881–1887, 1974. 
[68] Oseen. Arkiv matematik astron. Fysik A19, 1, 1925. 
[69] Frank. I. liquid crystals. on the theory of liquid crystals. Discuss. Faraday Soc., 25:19–28, 1958. doi: 10.1039/DF9582500019. URL http://dx.doi.org/10.1039/ <BR>DF9582500019. 
[70] Gramsbergen, Longa, and de Jeu. Landau theory of the nematic-isotropic phase transition. Physics Reports, 135(4):195–257, 1986. 
[71] LongaandPająk.Luckhurst–romanomodelofthermotropicbiaxialnematicphase. Liquid crystals, 32(11-12):1409–1417, 2005. 
[72] Allender and Longa. Landau–de gennes theory of biaxial nematics reexamined. Physical Review E, 78(1):011704, 2008. 
[73] Greiner, Neise, Stöcker, and Rischke. Thermodynamics and statistical mechanics, volume 1. Springer-Verlag New York, 1995. 
[74] Huang. Introduction to statistical physics. CRC Press, 2001. 
[75] MetropolisandUlam.Themontecarlomethod. Journal of the American statistical association, 44(247):335–341, 1949. 
[76] Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller. Equation of state calculations by fast computing machines. The journal of chemical physics, 21:1087, 1953. 
[77] Heard. Rigid body mechanics. Wiley-VCH, 2008. 
[78] Marsaglia. Choosing a point from the surface of a sphere. The Annals of Mathematical Statistics, 43(2):645–646, 1972. 
[79] Podlozhnyuk. Parallel mersenne twister. NVIDIA white paper, 2007. 
[80] Durst. Using linear congruential generators for parallel random number generation. In Proceedings of the 21st conference on Winter simulation, pages 462–466. ACM, 1989. 
[81] Hortensius,McLeod,andCard.Parallelrandomnumbergenerationforvlsisystems using cellular automata. Computers, IEEE Transactions on, 38(10):1466–1473, 1989. 
[82] Marsaglia. 	Yet another rng, 1994. URL https://groups.google.com/forum/ <BR>?fromgroups=#!topic/sci.stat.math/p7aLW3TsJys. 
[83] Goresky and Klapper. 	Efficient multiply-with-carry random number generators with maximal period. ACM Transactions on Modeling and Computer Simulation (TOMACS), 13(4):310–321, 2003. 
[84] Gratton. Parallel random number generation. URL http://www.ast.cam.ac.uk/ <BR>~stg20/cuda/random/. 
[85] Brink and Satchler. Angular momentum. Clarendon Press Oxford, 1993. 
[86] Jerphagnon, Chemla, and Bonneville. The description of the physical properties of condensed matter using irreducible tensors. Advances in Physics, 27(4):609–650, 1978. 
[87] Mulder. Isotropic-symmetry-breaking bifurcations in a class of liquid-crystal models. Physical Review A, 39(1):360, 1989. 
[88] Acharya, Primak, Dingemans, Samulski, and Kumar. 	The elusive thermotropic biaxial nematic phase in rigid bent-core molecules. Pramana, 61(2):231–237, 2003. 
[89] Olivares, Stojadinovic, Dingemans, Sprunt, and Jákli. Optical studies of the nematic phase of an oxazole-derived bent-core liquid crystal. Physical Review E, 68 (4):041704, 2003. 
[90] Shimbo, Gorecka, Pociecha, Araoka, Goto, Takanishi, Ishikawa, Mieczkowski, Go-mola, and Takezoe. Electric-field-induced polar biaxial order in a nontilted smectic phase of an asymmetric bent-core liquid crystal. Physical review letters, 97(11): 113901, 2006. 
[91] Francescangeli, Stanic, Torgova, Strigazzi, Scaramuzza, Ferrero, Dolbnya, Weiss, Berardi, Muccioli, et al. Ferroelectric response and induced biaxiality in the nematic phase of bent-core mesogens. Advanced Functional Materials, 19(16):2592– 2600, 2009. 
[92] Nagaraj, Panarin, Manna, Vij, Keith, and Tschierske. Electric field induced biaxiality and the electro-optic effect in a bent-core nematic liquid crystal. Applied Physics Letters, 96(1):011106–011106, 2010. 
[93] Dunmur, Szumilin, and Waterworth. Field-induced biaxiality in nematics. Molecular Crystals and Liquid Crystals, 149(1):385–392, 1987. 
[94] Barberi, Ciuchi, Durand, Iovane, Sikharulidze, Sonnet, and Virga. Electric field induced order reconstruction in a nematic cell. The European Physical Journal E, 13(1):61–71, 2004. 
[95] Borshch, Shiyanovskii, and Lavrentovich. Electric field induced biaxial order and differential quenching of uniaxial fluctuations in a nematic with negative dielectric anisotropy. Molecular Crystals and Liquid Crystals, 559(1):97–105, 2012. 
[96] Palffy-Muhoray and Dunmur. 	Field-induced order in nematic liquid crystals. Molecular Crystals and Liquid Crystals, 97(1):337–349, 1983. 
[97] Remler and Haymet. Phase transitions in nematic liquid crystals: a mean-field theory of the isotropic, uniaxial, and biaxial phases. The Journal of Physical Chemistry, 90(21):5426–5430, 1986. 
[98] Carlsson and Leslie. Behaviour of biaxial nematics in the presence of electric and magnetic fields: Evidence of bistability. Liquid Crystals, 10(3):325–340, 1991. 
[99] Berardi, Orlandi, and Zannoni. Columnar phases and field induced biaxiality of a gay–berne discotic liquid crystal. Physical Chemistry Chemical Physics, 2(13): 2933–2942, 2000. 
[100] Berardi, Muccioli, and Zannoni. Field response and switching times in biaxial nematics. The Journal of chemical physics, 128(2):024905–024905, 2008. 
[101] Peroukidis, Karahaliou, Vanakaras, and Photinos. Biaxial nematics: symmetries, order domains and field-induced phase transitions†. Liquid Crystals, 36(6-7):727– 737, 2009. 
[102] Jose, Jayasri, Dhara, Murthy, and Sastry. Effect of external field on the nematic to isotropic transition: An entropic sampling study. Molecular Crystals and Liquid Crystals, 545(1):168–1392, 2011. 
[103] Trojanowski, Allender, Longa, and Kuśmierz. Theory of phase transitions of a biaxial nematogen in an external field. Molecular Crystals and Liquid Crystals, 540(1):59–68, 2011. 
[104] Biscarini, Chiccoli, Pasini, Semeria, and Zannoni. Phase diagram and orientational order in a biaxial lattice model: A monte carlo study. Physical review letters, 75 (9):1803–1806, 1995. 
[105] Chiccoli, Pasini, Semeria, and Zannoni. A detailed monte carlo investigation of the tricritical region of a biaxial liquid crystal system. International Journal of Modern Physics C, 10(02n03):469–476, 1999. 
[106] Sonnet, Virga, and Durand. Dielectric shape dispersion and biaxial transitions in nematic liquid crystals. Physical Review E, 67(6):61701, 2003. 
[107] Bisi, Virga, Gartland Jr, De Matteis, Sonnet, and Durand. Universal mean-field phase diagram for biaxial nematics obtained from a minimax principle. Physical Review E, 73(5):051709, 2006. 
[108] Longa, Grzybowski, Romano, and Virga. Minimal coupling model of the biaxial nematic phase. Physical Review E, 71(5):051714, 2005. 
[109] Allender, Lee, and Hafiz. Landau theory of biaxial and uniaxial nematic liquid crystals. Molecular Crystals and Liquid Crystals, 124(1):45–52, 1985. 
[110] Meyer. 	Piezoelectric effects in liquid crystals. Physical Review Letters, 22(18): 918–921, 1969. 
[111] Hahn. 	Cuba—a library for multidimensional numerical integration. Computer Physics Communications, 168(2):78–95, 2005. 
[112] Swendsen. 	How the maximum step size in monte carlo simulations should be adjusted. Physics Procedia, 15:81–86, 2011. 
[113] Helfrich. Effect of electric fields on the temperature of phase transitions of liquid crystals. Physical Review Letters, 24(5):201–203, 1970. 
[114] Rosenblatt. Magnetic field dependence of the nematic-isotropic transition temperature. Physical Review A, 24:2236–2238, 1981. 
[115] Rosenblatt. Search for a field-induced nematic-isotropic critical point. 	Physical Review A, 25:1239–1242, 1982. 
[116] Nicastro and Keyes. 	Electric-field-induced critical phenomena at the nematicisotropic transition and the nematic-isotropic critical point. Physical Review A, 30(6):3156, 1984. 
[117] Walba, Körblova, Shao, Maclennan, Link, Glaser, and Clark. A ferroelectric liquid crystal conglomerate composed of racemic molecules. Science, 288(5474):2181– 2184, 2000. 
[118] LongaandPająk.Generalizeddispersionmodeloforientationallyorderedphasesof bent–core liquid crystals. Molecular Crystals and Liquid Crystals, 541(1):152–390, 2011. 
[119] Trojanowski, Pająk, Longa, and Wydro. Tetrahedratic mesophases, chiral order, and helical domains induced by quadrupolar and octupolar interactions. Physical Review E, 86(1):011704, 2012. 
[120] Brand and Pleiner. Macroscopic behavior of non-polar tetrahedratic nematic liquid crystals. The European Physical Journal E, 31(1):37–50, 2010. 
[121] Fel. Tetrahedral symmetry in nematic liquid crystals. Physical Review E, 52(1): 702, 1995. 
[122] Romano. Computer simulation study of a simple tetrahedratic mesogenic lattice model. Physical Review E, 77(2):021704, 2008. 
[123] Swendsen and Wang. Replica monte carlo simulation of spin-glasses. 	Physical Review Letters, 57(21):2607–2609, 1986.