Robert A. Skutnik, Immanuel S. Geier, Martin Schoen
A biaxial nematic liquid crystal composed of
matchbox-symmetric molecules
Open Access via institutional repository of Technische Universität Berlin
Document type
Journal article | Accepted version
(i. e. final author-created version that incorporates referee comments and is the version accepted for
publication; also known as: Author’s Accepted Manuscript (AAM), Final Draft, Postprint)
This version is available at
https://doi.org/10.14279/depositonce-14849
Citation details
Robert A. Skutnik, Immanuel S. Geier & Martin Schoen (2020) A biaxial nematic liquid crystal composed of
matchbox-symmetric molecules, Molecular Physics, 118:9-10,
https://doi.org/10.1080/00268976.2020.1726520.
This is an Accepted Manuscript of an article published by Taylor & Francis in Molecular Physics on 11 Feb
2020, available online: http://www.tandfonline.com/10.1080/00268976.2020.1726520.
Terms of use
This work is protected by copyright and/or related rights. You are free to use this work in any way permitted by
the copyright and related rights legislation that applies to your usage. For other uses, you must obtain
permission from the rights-holder(s).
Abiaxialnematicliquidcrystalcomposedofmatchbox-symmetric
molecules
Robert A. Skutnika,ImmanuelS.Geier
a,andMartinSchoen
a,b
aStranski-Laboratorium für Physikalische und TheoretischeChemie,TechnischeUniversität
Berlin, Straße des 17. Juni 115, 10623 Berlin, Germany; bDepartment of Chemical
Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ,
United Kingdom
ARTICLE HISTORY
Compiled January 30, 2020
ABSTRACT
By means of Monte Carlo simulations in the isothermal-isobaric ensemble we investi-
gate the structure and phase behaviour of a thermotropic liquid crystal composed of
matchbox-symmetric (or board-like) molecules. Besides theisotropicphasetheliquid
crystal exhibits also uniaxial and biaxial nematic phases. The interaction potential is
derived through an expansion in terms of Stone’s rotational invariants [A. J. Stone,
Mol. Phys. 78,241–256,(1978)]thatcanbereexpressedintermsofCartesianten-
sors. This latter formulation is particularly well suited for computer simulations. We
analyse the orientation distribution function which allowsustodistinguishbetween
intrinsic and extrinsic biaxiality. In addition, we study the orientation-dependent
correlation functions. In the limit of large intermolecularseparationsthevalueof
the orientation correlation function corresponds to the uniaxial and biaxial order
parameters which are coupled in a complex fashion.
KEYWORDS
Liquid crystal; nematic phase; biaxiality; Monte Carlo simulation
1. Introduction
Liquid crystals are a fascinating class of soft-matter systems capable of forming a host
of ordered phases [1]. They consist of organic molecules (i.e., mesogens) where shape
and flexibility of the molecular skeleton in conjunction withfavourablethermodynamic
conditions promote the formation of these ordered phases. The simplest one is the so-
called nematic phase in which the distribution of the centres-of-mass of the mesogens
is always isotropic. Yet, the mesogens align to a certain extent with a distinguished
direction, the so-called nematic director !
n.Thenematicphaseobviouslyhasuniaxial
symmetry in this case. However, this is not necessarily always so.
For example, in 1970 Freiser [2] and a little later Alben [3] hypothesised that it
might be possible for nematic phases to exhibit biaxial instead of uniaxial symmetry.
This conjecture rests upon a mean-field approach in both cases. In order for a nematic
phase to exhibit biaxial symmetry, a second symmetry axis besides !
nneeds to exist.
In the realm of modern display technology, biaxial nematic liquid-crystalline materials
Corresponding author: R. A. Skutnik, e-mail: robert.skutnik@campus.tu-berlin.de
are very interesting because of their short response times toappliedexternalfields[4].
Usually, one distinguishes between intrinsic and extrinsic biaxiality [5]. Biaxiality is
intrinsic if the mesogens themselves possess two symmetry axes. On the contrary, biax-
iality is extrinsic if uniaxially symmetric mesogens are manipulated by external agents
so that the system as a whole can exhibit a second symmetry axisunderfavourable
conditions [6]. The external agent could be a second uniaxialcompoundinabinary
mixture where the mesogens pertaining to different components prefer a T-shaped ar-
rangement [7]. The same philosophy has been chosen in the workbyCuetoset al.
[8]. However, unlike Ref. [7], where orientation-dependentLennard-Jones-typeofpo-
tentials were employed, Cuetos et al. base their study on (discontinuous) square-well
interactions.
The first experimental realisation of a (lyotropic) biaxial nematic liquid crystal was
reported by Yu and Saupe [9]. They studied a ternary mixture ofpotassium-laureate-
1-decanol and water. Another lyotropic liquid crystal that exhibits biaxial nematic
phases is a ternary mixture of soap, detergent, and water studied by Oliveira et al. [10].
Lyotropic liquid crystals exhibiting extrinsic biaxialityhavebeenstudiedbyStroobants
and Lekkerkerker [11] who investigated a mixture of rod- and platelet-like mesogens
within the framework of Onsager’s theory. Using again Onsager’s theory, Wensink et
al. [12] report a first-order phase transition during which a stable biaxial nematic forms
in a lyotropic liquid crystal.
As far as nematic phases with intrinsic biaxiality are concerned quite a bit of progress
has been made because of advances in chemical synthesis. There are essentially two
main molecular geometries that have been studied, namely bent-core [13, 14] and board-
or matchbox-symmetric mesogens [15, 16]. Both molecular architectures differ in the
number of reflection planes. Whereas there are two such planesforbent-corestructures,
matchbox-shaped mesogens possess three of those. Slightly more extravagant architec-
tures are those of organosiloxane tetrapodes [17] where matchbox-shaped central units
are connected to a silicon atom. Because of this structure evenbiaxialshort-rangeor-
der could be detected by means of deuterium nuclear magnetic resonance (NMR). In
afewcasesevenbiaxialsmecticphaseshavebeenreported[18, 19].
On the theoretical side quite a bit of work has already been devoted to biaxial ne-
matic phases in liquid crystals. Tjipto-Margo and Evans [20]employedOnsager’stheory
to investigate the formation of uniaxial nematic phases in a liquid crystal composed of
biaxial mesogens. At its core, Onsager’s theory requires a calculation of the excluded
volume of hard bodies. This problem has been investigated in depth by Mulder [21].
By means of computer simulations Allen [22] and later Camp and Allen [23] inves-
tigated the phase behaviour of intrinsically biaxial mesogens with only hard repulsive
interactions. They proposed a set of order parameters also calculated in this work.
These order parameters are the same as those employed also in the work by Mulder
[21] and Straley [24].
Berardi et al. [25, 26] use a properly modified Gay-Berne potential in their sim-
ulations of a liquid crystal with intrinsic biaxiality. In a lattice Monte Carlo (MC)
study Biscarini et al. [27] used an orientation dependence of the interaction poten-
tial based upon symmetry-adapted Wigner rotation matrices.Theirmodelreducesto
the Lebwohl-Lasher model in the limit of vanishing biaxiality. Results from mean-field
theory are in good qualitative agreement with their MC data.
Recently, Cuetos et a. [28] employed MC simulations of board-like particles. These
particles interact through purely entropic forces. The authors focus on the impact of the
variation of the width-to-length ratio while keeping the thickness of the mesogens fixed.
They found a rich phase behaviour including biaxial columnarandsmecticphases.They
2
also use Onsager’s theory to predict the uniaxial-biaxial nematic phase transition.
Alatticemodelinvokinganexpansionoftheinteractionpotential based upon rota-
tional invariants [29] has been proposed by Luckhurst and Romano [30]. Their model
potential was used later by Longa and Pająk [31] who employed the Ginzburg-Landau
phenomenological theory to study the phase behaviour of mesogens with intrinsic bi-
axiality.
Amodelformatchbox-symmetricmesogenswasusedwithinadensity functional ap-
proach by Longa et al. [32]. The impressive richness of the phase diagram of matchbox-
symmetric hard mesogens was demonstrated by Martínez-Ratón[33].Theirstudyis
based upon fundamental measure theory for the contribution of the hard-repulsive
interactions to the excess free energy.
In this work we utilise an off-lattice model for a liquid crystal with intrinsic biax-
iality. The molecular architecture assumes matchbox-symmetric mesogens. To obtain
an expression for the orientation dependence of the anisotropic interactions we fol-
low Luckhurst and Romano [30] and expand the interaction potential in terms of the
rotational invariants proposed by Stone [34].
Based upon symmetry considerations [21, 34] and utilising a ladder-operator for-
malism [29] borrowed from angular-momentum theory in quantum mechanics [35], we
can derive an expression for the anisotropic interaction potential in terms of a full
contraction of Cartesian second-rank tensors. Thus, our potential is equivalent to the
one used in the mean-field calculations of Straley [24] and of Sonnet et al. [5] (see also
Ref. 36). We focus on uniaxial and biaxial order parameters [22, 23], the orientation
distribution function (odf), and the orientation correlation function.
The odf reflects the symmetry of the mesogens in our model. In addition, we show
that at large separations between the mesogens the orientation correlations can be
expressed in terms of the uniaxial and biaxial order parameters with a subtle coupling
between the two. Such a coupling obviously does not occur in uniaxial nematic phases
where much simpler expressions are obtained [37]. Neither the odf nor the orientation
correlation function have been considered in earlier work for similar systems. Both
quantities turn out to provide a lot of insight into phase changes in this model.
The remainder of the manuscript is organised as follows. Our model is introduced
in Section 2. Important mathematical details of the derivation of the anisotropic con-
tribution to the interaction potential are deferred to Appendix A. We give details of
the mathematical manipulations here because they are difficult to extract from the
existing literature. In Section 3 we introduce key quantities that we intend to analyse
in the subsequent Section 4. The paper concludes in Section 5 where we put our re-
sults into perspective. Details of the analysis of the asymptotic decay of the orientation
correlation functions are given in Appendix B.
2. Model
We consider a liquid crystal composed of Nmesogens. The interaction between meso-
gens 1and 2,whosecentresofmassarelocatedatr1and r2,isassumedtobepairwise
additive. We decompose the interaction potential uinto an isotropic part uiso and into
an anisotropic contribution uaniso according to
u(r12,Ω1,Ω2)=uiso(r12)+uaniso(r12,Ω1,Ω2),(1)
3
where r12 =|r12|=|r1−r2|is the distance between the centres of mass. In Equa-
tion (1), Ω1,2are Euler angles. We adopt the convention of Gray and Gubbins [38]
for the order in which the rotations about specific axes are being carried out (see Sec-
tion A.2 of Ref 38). Therefore, Ω=(φθχ)where φand θrefer to the azimuthal and
polar angle, respectively.
For uiso we adopt the soft-sphere potential which permits us to express uiso as
uiso(r12)=4ε"σ
r12 #12
,(2)
where, for the time being, εand σare just two parameters that will be given a lucid
interpretation shortly. To derive an expression for uaniso we follow Stone [34] who
noted that for molecules of arbitrary shape an expansion of uaniso in a complete set of
rotational invariants according to
uaniso (r12,Ω1,Ω2)=$
l1l2l$
n1n2
un1n2
l1l2l(r12)
×Sn1n2
l1l2l(Ω1Ω2,Ω12),(3)
is generally feasible. In Equation (3), {un1n2
l1l2l}is a set of expansion coefficients that
depend only on r12,andΩ12 =(φ12θ120) describes the orientation of !
r12 =r12/r12
in the space-fixed reference frame where !
r12 is pointing from the centre of mass of
mesogen 1to that of mesogen 2.Hereandbelowweusethecarettoindicateaunit
vector. The quantity l!(that is l1,l2,orl)isanon-negativeinteger(zeroincluded)
and m!and n!are related to l!through the relation m!,n
!∈[−l!,l!].
ui
viwi
v
wu
h
Figure 1. Cartoon of a matchbox-symmetric mesogen for which the major symmetry axis !
uiand the minor
symmetry axis !
viare given in Equations (11a) and (11b), respectively; the direction of the second minor
symmetry axis is obtained from !
wi=!
ui×!
vi.Themesogencontainsreflectionplanesσhand σwu
vin the
notation of Stone (see Table 2 of Ref. 34). Both reflection planes contain the vector !
wi.Athirdreflectionplane
σvu
v,orthogonaltoσhand σwu
v,isnotshown.
4
Table 1. Val-
ues of l,n1,and
n2to be con-
sidered in Equa-
tion (6).
ln
1n2
00 0
20 0
±20
0±2
±2±2
Moreover, in Equation (3) we introduce the scalar function [34]
Sn1n2
l1l2l(Ω1,Ω2,Ω12)=il1−l2−l$
m1m2m"l1l2l
m1m2m#
×Dl1
m1n1
∗(Ω1)Dl2
m2n2
∗(Ω2)Dl
m0
∗(Ω12),(4)
where i≡√−1,(···)is a Wigner 3jsymbol,
Dl
mn(Ω) = e−imφdl
mn (θ)e−inχ(5)
is a Wigner rotation matrix, and the asterisk denotes the complex conjugate [35, 38, 39].
The real matrix dl
mn depends only on the Euler angle θ;forl≤2it is given by
Equations (A.114) and (A.115) of the book by Gray and Gubbins [38]. The asymmetric
phase factor in front of the summation signs serves to make Sn1n2
l1l2lreal for n1=n2=0
regardless of l1,l2,andl[34].
We now specialise the above treatment to matchbox-symmetricparticlesillustrated
by the cartoon presented in Figure 1. Each mesogen possesses three symmetry axes:
the major one labelled !
uiand two minor ones that we refer to as !
viand !
wi.These
three axes are obviously pairwise orthogonal to one another. Because of the matchbox
symmetry, each mesogen contains three reflection planes. Twooftheselabelledσhand
σwu
vin Figure 1 contain the minor symmetry axis !
wi.AccordingtoFigure1,!
uiis a
C2rotation axis (and so are !
viand !
wi).Therefore,thematchbox-symmetricmesogens
pertain to the point group D2h.
The existence of these symmetry elements has consequences for the expansion coeffi-
cients in Equation (3). From Figure 1 one readily sees that each mesogen is centrosym-
metric. According to Stone [34], lin Equation (4) must therefore be even. Moreover,
the existence of the reflection plane σh(see Figure 1) then requires n1and n2to be
even as well [34].
In addition, we assume uaniso to depend only on r12.ThisimpliesthatinEqua-
tion (4), l=m=0. Because of this assumption, the 3jsymbol in Equation (42) is
nonzero only for l1=l2and m1=−m2.Thisfollowsfromthetwoselectionrules
•|l1−l2|≤l≤l1+l2(triangle inequality) and
•m1+m2=m
which need to be satisfied simultaneously. On account of theseconsiderationsthetriple
sum in Equation (4) collapses to a single one; the same is true for the triple sum in
Equation (3).
However, the sum on min Equation (4) involves 2l+1 summands and thus becomes
5
overwhelmingly large rather quickly as lincreases. Hence, we restrict the discussion to
l≤2from now on and rewrite Equation (3) as
uaniso (r12,Ω1,Ω2)=
2
$
l=0 $
n1n2
un1n2
ll0(r12)Sn1n2
ll0(Ω1,Ω2,Ω12)(6)
keeping the argument Ω12 of the %Sn1n2
ll0&only formally. Hence, in principle we need to
consider in Equation (6) the terms corresponding to the entries in Table 1. However,
the ten terms listed in Table 1 can be reduced by additional symmetry considerations
[21, 30] as we explain in the following.
We begin with the simplest case l=n1=n2=0. Because S00
000 =1,theleading
term in the expansion in Equation (6) is u00
000.However,wecannottakerecoursetoany
additional first-principles argument to determine what thisexpansioncoefficientshould
be. We therefore assume that our mesogens interact via dispersion forces. Under this
premise, Gray and Gubbins [38] pointed out that u00
000 ∝r−6
12 .Thecoefficientofpro-
portionality has dimensions of energy ×(length)−6[see, for example, Equation (2.226)
of Ref. 38]. This prompts us to introduce
u00
000(r12)=−4ε"σ
r12 #6
.(7)
Moreover, from Equation (2.234) of Ref. 38 it is clear that alltheotherexpansion
coefficients of uaniso should have the same distance dependence as long as we restrict
the discussion to non-polarisable mesogens interacting viadispersionforcesonly.
Next, we consider l=2but n1=n2=0.AsweshalldemonstrateinAppendixA
this contribution refers to a coupling between purely uniaxial degrees of freedom which
we shall keep. The next two cases, namely l=2,n1=±2and n2=0as well as
n1=0and n2=±2are disregarded because they would entail a coupling between
uniaxial and biaxial degrees of freedom as we shall demonstrate also in Appendix A.
We deliberately omit such coupling in the present work.
Finally, we include the cases l=2,n1=±2,n2=±2.Thenotationindicates
that all four combinations of negative and positive values ofn1and n2have to be
considered. The corresponding terms in our expression for uaniso can be simplified by
additional symmetry considerations. For example, according to Table 2 in the paper
by Stone [29],
u22
220 (r12)=u22
220 (r12)(8a)
u22
220 (r12)=u22
220 (r12),(8b)
where we used the shorthand notation x=−x.TheseparateequalitiesgiveninEqua-
tions (8a) and (8b) follow from the existence of the reflectionplaneσwu
v(see Figure 1)
[21, 30]. Reminding ourselves that the mesogens are identical, the expansion coefficients
in Equation (6) satisfy one last equality namely u22
220 =u22
220 from which it follows that
all four expansion coefficients in the Equations (8) are equal. Because of these addi-
tional symmetry considerations the number of entries listedinTable1isfinallyreduced
to only 5.Thus,forthecasen1=±2,n2=±2the expansion coefficients in Equa-
tion (6) are all the same. This allows us to combine the associated rotational invariants
6
and introduce
S±2±2
220 =S22
220 +S22
220 +S22
220 +S22
220.(9)
Rather than expressing uaniso in terms of rotational invariants {Sn1n2
220 }it will turn out
to be more useful to express the anisotropic interaction potential in terms of suitably
defined Cartesian tensors. Let us now introduce the traceless1,second-ranktensors
U(Ωi)≡Ui=!
u(Ωi)!
u(Ωi)−1
31(10a)
B(Ωi)≡Bi=!
v(Ωi)!
v(Ωi)−!
w(Ωi)!
w(Ωi),i=1,2,(10b)
where 1is the unit tensor and
!
u(Ωi)≡!
ui=
sin θicos φi
sin θisin φi
cos θi
(11a)
!
v(Ωi)≡!
vi=
cos φicos θicos χi−sin φisin χi
sin φicos θicos χi+cosφisin χi
−sin θicos χi
(11b)
represent the orientation of the molecular major and minor symmetry axes in a space-
fixed coordinate system. The tensors Uand Bare referred to as the uniaxial and biaxial
tensors, respectively. The terms involving products of !
u,!
v,and !
won the right-hand
sides of Equations (10a) and (10b) are to by understood as dyads.
The third vector is defined through the expression !
wi=!
ui×!
vi.Usingtheexpressions
in Equations (2) and (6)–(10) we can eventually cast uin Equation (1) as
u(r12,Ω1,Ω2)=uLJ (r12)−4ε"σ
r12 #6
Ψ(Ω
12),(12)
where
uLJ =4ε+"σ
r12 #12
−"σ
r12 #6,(13)
is the well-known Lennard-Jones potential,
Ψ(Ω
12)=ε!U1:U2+ε!!B1:B2(14)
is the anisotropy function, and X:X=XαβXβα (X=Uor B)denotesthefull
contraction of two second-rank tensors using the Einstein repeated summation-index
convention.
At this point we emphasise that the expressions given in Equations (12)–(14) do by
no means allow for modeling mesogens that have the shape of matchboxes. Rather the
interaction potential possesses certain symmetry elements such as n-fold rotation axes
and reflection planes that would be characteristic of a real matchbox. Therefore, we
refer to our mesogens as “matchbox-symmetric”.
1For particularly important properties of traceless matrices, see Ref. [40].
7
In Appendix A we link the rotational invariants introduced inEquation(4)tothe
Cartesian tensors Uand Bin Equations (10a) and (10b), respectively. This is accom-
plished with the ladder-operator technique proposed by Stone [34]. It is identical to
the one arising in angular-momentum theory in quantum mechanics [35].
The philosophy of just equipping the interaction potential with certain symmetry
elements is not at all new. In fact, it has been employed in the original work by
Maier and Saupe [41] in which the orientation dependence of the interaction between
uniaxial mesogens is accounted for by the second Legendre polynomial. As pointed out
by Schoen et al. [42], for any fixed relative orientation of a pair of mesogens,Maier
and Saupe’s interaction potential will be isotropic. Nonetheless, properties of uniaxial
nematic phases have been successfully reproduced by this minimalistic approach.
Finally, we are in a position to attach physical meaning to theparametersεand σ
introduced in Equation (2) through the relations uLJ (σ)=0and duLJ/dr12|r12 =σ=
0.Thus,σis a measure of the diameter of the spherically symmetric coreofour
mesogens and εis the depth of the attractive well. In Equation (14), ε!and ε!! are
dimensionless coupling constants referring to the strengthofthecouplingofuniaxial
and biaxial degrees of freedom, respectively. In these parameters we lump together
all trivial numerical pre-factors arising from the expansion of uaniso in terms of the
rotational invariants %Sn1n2
l1l2l&and other molecular constants such as polarizabilities
[38].
Throughout this work we fix ε!=0.375.Thisvaluewaschosenaccordingtothe
following criteria. On the one hand, ε!must not be too large to avoid glassy states
before biaxial nematic phases can even form. On the other hand, ε!must not be too
small so that stable uniaxial and biaxial nematic phases occur for moderate values of T
and ρ.Inpractice,weperformedshorttestrunstofine-tunethisparameter according
to both criteria.
3. Orientational structure and order parameters
3.1. Orientation-dependent correlation functions
3.1.1. Expansion coefficients
To gain deeper insight into the nature of uniaxial and biaxialnematicorderweintro-
duce the orientation-dependent pair correlation function.Asourpointofdeparturewe
choose the two-particle generic distribution function [38]definedas
ρ(r1,r2,Ω1,Ω2)=ρ(r1,Ω1)ρ(r2,Ω2)g(r1,r2,Ω1,Ω2),(15)
where gis the orientation-dependent pair correlation function andρon the right hand-
side denote the generic single-particle distribution functions. Assuming our system to
be homogeneous we can rewrite the latter function as
ρ(r,Ω) = ρα (Ω),(16)
where ρon the right-hand side is the (mean) number density and αis the odf. The
odf is normalised according to -dΩα=1. Because homogeneity implies translational
8
invariance we can rewrite Equation (15) as
1
ρ2ρ(r12,Ω1,Ω2)=α(Ω1)α(Ω2)g(r12,Ω1,Ω2)
=g!(r12,Ω1,Ω2),(17)
where g!is a “reduced” distribution function [37]. It is introduced because in the absence
of positional and orientational correlations gr12→∞
−→ 1and therefore g!can simply be
cast as the product of the two odf’s [37].
Similar to Equation (3) we can expand g!in terms of rotational invariants via
g!(r12,Ω1,Ω2)=$
l1l2l$
n1n2
gn1n2
l1l2l(r12)Sn1n2
l1l2l(Ω1,Ω2,Ω12),(18)
where {gn1n2
l1l2l}is a set of expansion coefficients. We immediately restrict thediscussion
to the case in which gdepends only on r12.Thisimpliesthatl=m=0.Employing
the two selection rules for nonzero 3jsymbols introduced in Section 2 we can replace
the 3jsymbol by [see Equation (4)]
"ll0
mm0#=(−)l+m
√2l+1.(19)
Using Equation (4) we can rewrite Equation (18) more explicitly as
g!(r12,Ω1,Ω2)=$
l
1
√2l+1"2l+1
8π2#2$
n1n2
gn1n2
ll0(r12)
×$
m
(−)mDl
mn1
∗(Ω1)Dl
mn2
∗(Ω2),(20)
where (−)m=(−1)m.
The extra factor of .2l+1
8π2/2in Equation (20) has been introduced to make sure
that the final expression for gn1n2
ll0is compatible with the Cartesian tensor formulation
developed in Appendix A. The latter does not take recourse to the orthogonality of
the Wigner rotation matrices and thus lacks factors of this sort.
Multiplying both sides of Equation (20) by
$
m
(−)mDl
mn1(Ω1)Dl
mn2(Ω2),
integrating over dΩ1dΩ2,andemployingtheorthogonalityoftheWignerrotationma-
trices
0dΩDl
mn
∗(Ω)Dl!
m!n!(Ω) = 8π2
2l!+1δl!lδm!mδn!n,(21)
where the δij’s are Kronecker symbols.
9
After rearranging terms we then obtain the relation
√2l+1gn1n2
ll0(r12)=00dΩ1dΩ2g!(r12,Ω1,Ω2)
×$
m
(−)mDl
mn1(Ω1)Dl
mn2(Ω2).(22)
Using Equation (4) it is easy to verify that our final expression for the expansion
coefficients gn1n2
ll0can be cast as
gn1n2
ll0(r12)=00dΩ1dΩ2g!(r12,Ω1,Ω2)Sn1n2
ll0
∗.(23)
3.1.2. Computer simulations
Following Zannoni [37, 43], we can rewrite Equation (23) to make it suitable for an
evaluation in computer simulations. To that end we divide andmultiplyby
g!
0(r12)=00dΩ1dΩ2g!(r12,Ω1,Ω2)(24)
which permits us to recast Equation (23) as [42, 44]
gn1n2
ll0(r12)=g!
0(r12)1Sn1n2
ll0
∗2shell ,(25)
where )...*shell ≡g!
0
−1--dΩ1dΩ2g!... and we dropped the arguments of Sn1n2
ll0to
simplify the notation.
Notice that on account of the definition of g!in Equation (17), g!
0is not the normal
radial distribution function. For one, g!still involves the product of the two odf’s
and second, the expression in Equation (24) is not defined as an unweighted average
over orientations of gintroduced in Equation (15). Therefore, g!
0is unknown at this
stage because the odf’s are undetermined. However, in practice this is not a problem
because orientation correlations are fully accounted for by1Sn1n2
ll02shell.Thus,ifthe
focus is exclusively on orientational correlations we can consider 1Sn1n2
ll02shell rather
than the expansion coefficient gn1n2
ll0associated with it.
We can now invoke the same symmetry considerations already applied to {un1n2
ll0}in
Section 2. Leaving out the case l=n1=n2=0for obvious reasons, only l=2and
n1=n2=0and n1=±2and n2=±2will be considered (see Table 1); the other two
cases, namely n1=±2,n2=0as well as n1=0and n2=±2are disregarded because
Ψin Equation (14) does not contain a term corresponding to coupling between uniaxial
and biaxial degrees of freedom [see also Equations (A13) and (A14)]. Therefore, we only
consider [see Appendix A]
√51S00
220 (Ω1,Ω2,Ω12)2shell =3
2)U1:U2*shell (26a)
√51S±2±2
220 (Ω1,Ω2,Ω12)2shell =)B1:B2*shell ,(26b)
where the tensors Uand Bhave been defined in Equations (10a) and (10b), respec-
tively.
10
Acoupleofcommentsapplyatthispoint.Forn1=n2=0,S00
ll0is real and therefore
S00
ll0
∗=S00
ll0so that Equation (26a) follows from the general expression given in Equa-
tion (23) without further ado. Using the symmetry property oftheWignerrotation
matrices
Dl
mn(Ω) = (−)m+nDl
mn
∗(Ω)(27)
it follows that
$
m
(−)mDl
mn1(Ω1)Dl
mn2(Ω2)
=$
m
(−)mDl
mn1
∗(Ω1)Dl
mn2
∗(Ω2)
=$
m
(−)mDl
mn1
∗(Ω1)Dl
mn2
∗(Ω2),(28)
where we omitted a factor (−)2m+n1+n2because n1and n2are supposed to be even
integers. Thus, it is clear that Sn1n2
ll0
∗=Sn1n2
ll0and therefore S±2±2
220
∗=S±2±2
220 using
Equation (9).
Utilising Equations (26) has the additional advantage that in a computer simulation
it is rather straightforward to compute )Sn1n2
220 *shell.Noticethattheexpressionsin
Equations (26a) and (26b) clearly indicate that only the relative orientations between
mesogens 1and 2matter. As explained by Zannoni [43] one could therefore compute
)Sn1n2
220 *shell as a histogram. According to the distance r12 ahistogrambinischosen
and the tensors Uand Bare computed for both mesogens of the pair. These tensors
are then fully contracted and the resulting values are sortedintothebincorresponding
to r12.Theprocedureisrepeatedforeachpairofmesogensinthegiven configuration.
Obviously, this procedure has to be repeated for as many configurations as it may take
to compute the histogram with an acceptable statistical accuracy. In closing we note
that a procedure very similar to the one described has been employed by Streett and
Tildesley [44] and later by Schoen et al. [42].
3.1.3. The limit of large intermolecular separations
It is instructive to investigate the expressions given in Equations (26a) and (26b) in
the limit r12 →∞.Inthislimit,wecanrewriteEquation(24)as
lim
r12 →∞ g!
0(r12)=00dΩ1dΩ2α(Ω1)α(Ω2)=1,(29)
where we also used Equation (17). The previous result followsbecausetheodf’sare
normalised and because we assumed that limr12 →∞ g=1[37]. This hypothesis is
certainly valid in an isotropic phase; in an ordered phase, however, the assumption can
only be justified aposteriori.Thisisbecauselimr12 →∞ g=1neglects all (spatial and
orientational) correlations between mesogens 1and 2which is difficult to justify at the
outset.
However, accepting the above approximation as a working hypothesis for the time
11
being we use Equation (23) and introduce
lim
r12 →∞ gn1n2
220 (r12)=)Sn1n2
220 *∞
=00dΩ1dΩ2α(Ω1)α(Ω2)Sn1n2
220
∗
=1
√500dΩ1dΩ2α(Ω1)α(Ω2)
×$
m
(−)mD2
mn1(Ω1)D2
mn2(Ω2).(30)
We dropped the word “shell” as an index because in the limit r12 →∞the shell loses
its meaning. In fact, gn1n2
220 will assume a constant value irrespective of how one defines
ashell.
We now follow Mulder [21] and expand the odf according to
α(Ω) = $
l!m!n!
αl!
m!n!Ql!
m!n!(Ω)
=$
lmn
αl
|m||n|s|m||n|Dl
mn(Ω),(31)
where m!,n
!∈[0,l];instead,nand mare defined as usual (see Section 2). In addition,
the symmetry of the odf requires l,m,andnto be even. The functions %Ql
m!n!&are
defined in Equation (3.5) of the paper by Mulder [21]. The coefficient
smn ="1
√2#2+δm0+δn0
2δm0+δn0=2
(δm0+δn0−2)/2(32)
serves to make the expression on the second line of Equation (31) consistent with
Mulder’s Equations (3.5) and (3.6). The coefficients %αl
mn&are expansion coefficients
[21–24, 37, 43]. The expansion coefficients should not be confused with order parameters
because the latter are defined on the interval [0,1] whereas this is not necessarily true
for the former,
We now insert Equations (31) and (32) into Equation (30) and obtain eventually
[cf., Equations (26)]
√51S00
2202∞=3
2)U1:U2*∞=.α2
00/2+.α2
20/2(33a)
√51S±2±2
220 2∞=)B1:B2*∞=23.α2
2222+1α2
02/24,(33b)
where we dropped all arguments to simplify the notation. Details of this derivation
have been deferred to Appendix B. Our analysis also enables us to evaluate 1S±20
220 2∞
and 1S0±2
220 2∞where a definition analogous to the one for S±2±2
220 in Equation (9) is
used. Based upon Equations (A13) and (A14) we find that
√51S±20
220 2∞=53
2)U1:B2*∞=53
2)B1:U2*∞
=√26α2
00α2
02 +α2
22α2
207.(34)
12
Together the expressions in Equations (33) and (34) reflect the coupling between uniax-
ial and biaxial degrees of freedom. However, from the definition of Ψin Equation (14)
it is clear that we deliberately disregard such coupling. According to one’s physical
intuition one would therefore expect B1:U2=U1:B2=0even in uniaxial or bi-
axial nematic phases. If this is indeed so α2
20 =α2
02 =0and therefore we obtain from
Equations (33) our final result
√51S00
2202∞=S2(35a)
√51S±2±2
220 2∞=2η2,(35b)
where the nematic and biaxial order parameters Sand ηwill be defined below in
Equations (42) and (50). Equations (35) are derived by comparing Equation (3.9) in
Mulder’s paper with our Equations (42) and (50) and by noticing that our expansion
coefficients %αl
mn&have to be identical to Mulder’s expansion coefficients %ql
mn&[21].
We defer a study of the effects of including the coupling between uniaxial and biaxial
degrees of freedom to a separate publication.
3.2. Rotations of the eigensystems
Another quantity of interest to us is the odf. In computer simulations the odf can
be computed as a two-dimensional histogram of the angles θiand φidescribing the
orientation !
uiof mesogen i[see Equation (11a)]. Hence, the tips of the vectors {!
ui}
correspond to points on S2(i.e., the surface of the unit sphere). To compute the odf,
S2is partitioned into small bins of size δA=δφ ×δθ.Onethensimplycountsthe
number of mesogens whose orientation pertains to a specific bin. Clearly, the size of
the bins is largest along the equator of S2and becomes progressively smaller towards
its poles if the number of bins stays fixed.
Suppose now the liquid crystal is in a uniaxial nematic phase.Thenthevectors{!
ui}
are predominantly aligned with !
nwhere the tip of !
ncan be represented by a point
on S2as well. In the absence of external fields, this point can lie anywhere on S2.In
particular, it can be located near the poles. This implies that most of the {!
ui}also
correspond to points in that region of S2.Thus,ifthehistogramiscomputed,theodf
would only be accessible with rather low resolution, or, if the resolution is enhanced
by making δAsmaller, with insufficient statistical accuracy.
It is therefore advantageous to rotate the vectors {!
ui}(or, alternatively, the basis
used to define S2)suchthattheorientationsofthemesogensaresortedintobins
located at or near the equator of S2.Wecandetermineabasisoptimisedaccordingto
the above criteria following the protocol originally proposed by Camp and Allen [23]
to define suitable order parameters.
Let us introduce three instantaneous,second-rank,symmetric,andtracelesstensors
Qxx =1
2N
N
$
i=1
(3!
xi!
xi−1),(36)
where !
xirepresents !
ui,!
vi,or !
wi.ThethreetensorsQxx satisfy the eigenvalue equation
Qxx !
qx
±,0=λx
±,0!
qx
±,0,(37)
where λx
−≤λx
0<λ
x
+are the three eigenvalues and !
qx
±,0are the associated eigenvectors
13
where the equal sign applies if the nematic phase is uniaxial.Equation(37)givesus
atotalofnineeigenvaluesandassociated-vectors.Ithasbeen pointed out by Mulder
[21] that on account of symmetry arguments the three tensors Quu,Qvv and Qww share
their eigensystem and can thus be diagonalised simultaneously. However, the role each
of the three eigenvectors is playing depends on the symmetry axis one is referring to.
This is illustrated by the sketch in Figure 2 where the polar angle θvis formed between
!
viand the axis !b3;theazimuthalangleϕvis formed between the projection of !
vionto
the !b1–!b2plane and the !b1-axis.
b2
b1
,
b1b3
,
b3b2
,
ui'
vi'
v
v
uu
Figure 2. Sketch of the eigensystem (!
b1,!
b2,!
b3)(see text). Blue symbols are defined with respect to !
v!
i
whereas symbols in red refer to !
u!
i.Polarandazimuthalanglesθvand ϕvof the symmetry axis !
v!
iare defined
as indicated in the figure. The role played by the three eigenvectors changes if the symmetry axis !
u!
iis considered
as indicated such that !
b1→!
b3,!
b2→!
b1,and!
b3→!
b2(red symbols) and the corresponding angles θuand ϕu
are defined accordingly. The transformation (!
b1,!
b2,!
b3)→(!
b2,!
b3,!
b1)preserves the right-handedness of the
eigensystem.
We define as !b3the eigenvector corresponding to the eigenvalue λu
+.Similarly,we
take !b1to be the eigenvector corresponding to the eigenvalue λv
+.Lastbutnotleast,
!b2=!b3×!b1.Letusnowintroducethestandardbasis(!
ex,!
ey,!
ez),where!
eT
x=(1,0,0),
!
eT
y=(0,1,0),!
eT
z=(0,0,1) and Tdenotes the transpose. We can then transform the
vectors !
uiand !
vigiven in the standard basis to new vectors !
u!
iand !
v!
iusing the three
eigenvectors !b1,!b2,and!b3as the new reference frame. The transformations can be
cast as
!
u!
i=Ruu !
ui(38a)
!
v!
i=Rvv !
vi,i=1,...,N, (38b)
where !
u!
iand !
v!
idepend on ϕu,v and θu,v as illustrated by the sketch in Figure 2. Note,
that ϕu,v and θu,v themselves are functions of the three Euler angles φi,θi,andχi.
In Equations (38), Ruu and Rvv are respective direction-cosine matrices between the
vectors of the standard basis and (!b2,!b3,!b1)(reference to !
ui)ontheonehandand
(!b1,!b2,!b3)(reference to !
vi)on the other hand.
14
The relevant angles can then be computed via (see Figure 2)
cos θu,i =!
u!
i·!b1(39a)
cos θv,i =!
v!
i·!b3(39b)
and
cos ϕu,i =(!b2!b2+!b3!b3)!
u!
i·!b1
888(!b2!b2+!b3!b3)!
u!
i888
(40a)
cos ϕv,i =(!b1!b1+!b2!b2)!
v!
i·!b2
888(!b1!b1+!b2!b2)!
v!
i888
.(40b)
The tensors in parentheses in Equations (40a) and (40b) effectaprojectionof!
u!
iand
!
v!
ionto the !b2–!b3and !b1–!b2planes, respectively.
With the aid of Equations (39a) and (40a) on the one hand and (39b) and (40b) on
the other hand we are in a position to compute the odf from the expression
αu,v (ω)=
N
$
i=1 1δ.ω−ω!
i/2
2π
00
π
00
sin θdθdϕ
N
$
i=1 1δ.ω−ω!
i/2
,(41)
where )...*denotes an ensemble average [7], δis the Dirac δ-function, and δ(ω−ω!
i)is
shorthand notation for δ(θ−θ!
i)δ(ϕ−ϕ!
i)/sin θ.Thenotationαu,v is used to indicate
that we have to consider two odf’s depending on whether ω!
i≡ω!
u,i =(ϕu,i,θ
u,i)or
ω!
i≡ω!
v,i =(ϕv,i,θ
v,i)(see Figure 2).
There is an additional benefit of using the rotated coordinatesystems.Aswewill
demonstrate below the structure of the odf’s turns out to be particularly simple in the
rotated coordinate systems. For example, in a uniaxial nematic phase, αuwill exhibit
two isolated maxima because of the reflection plane σh(see Figure 1); at the same time,
αvwill have a band-like structure. In a biaxial nematic phase this band will change to
asequenceofisolatedmaximaseparatedbyacertainangleincrement.
3.3. Order parameters
To classify various disordered and ordered phases in our system we now introduce order
parameters via Equations (36) and (37). From the discussion in the preceding section
we saw that the (instantaneous) eigensystem formed by the three vectors !b1,!b2,and
!b3is shared by the three (instantaneous) tensors Quu,Qvv,andQww.Moreover,by
definition, !b1,!b2,and!b3define a rotating reference frame. This implies that during the
course of a simulation none of these rotating frames is favoured over any other including
that represented by the standard basis. This allows us to simplify the subsequent
analysis by assuming that !b1=!
ex,!b2=!
ey,and!b3=!
ezwithout loss of complete
generality.
15
We begin with the nematic order parameter defined as [cf., Ref.45]
S=1
20dΩ{3[!
u(Ω)·!b3]2−1}α(Ω).(42)
The integrand in this equation contains the scalar product (!
u·!b3)2squared which
reflects the existence of the reflection plane σh(see Figure 1). We can rewrite this
expression using the tensor identity
(!
u·!b3)2=[(
!
u·!b3)!
u]·!b3=!b3·!
u!
u·!b3.(43)
Noting that 1=!b3·1·!b3,wecanrewriteEquation(42)as
S=0dΩ[!b3·Q!
uu(Ω)·!b3]α(Ω),(44)
where [cf., Equation (36)]
Q!
uu(Ω) = 1
2[3!
u(Ω)!
u(Ω)−1].(45)
Alternatively, we may again start from Equation (42) and notethatitmaybe
rewritten as
S=0dΩD2
00(Ω)α(Ω),(46)
where we have used the fact that on account of our choice !b3=!
ezthe integrand can
be cast as
P2(cos θ)=1
2(3 cos2θ−1) = D2
00(Ω),(47)
where cos θ=!
u·!b3.UsinginEquation(46)theexpansionoftheodfintroducedin
Equation (31) we obtain
S=$
lmn
αl
|m||n|s!
|m||n|0dΩD2
00
∗(Ω)Dl
mn(Ω),(48)
where s!
mn =smn(2l+1)/8π2and smn is given in Equation (32). Invoking Equation (21)
it is straightforward to verify that Equation (48) can be recast as
S=α2
00.(49)
It should be noted though that we included an extra factor of (2l+1)/8π2in the
definition of s!
mn to ensure that the expansion coefficients {αl
|m||n|}become order pa-
rameters in the sense that the members of the set {αl
|m||n|}take values in the interval
[0,1].Wealsonoteinpassingthattogethertheexpressionsgivenin Equations (46)
and (47) agree with the one given by Camp and Allen [23] in theirEquation(17),
Equation (3.9) of Mulder [21], and Equation (5) in the paper byStraley[24].
16
Similarly, we can define the biaxiality order parameter through the expression
η≡1
20dΩ[(!b1·!
v)2+(
!b2·!
w)2
−(!b2·!
v)2−(!b1·!
w)2]α(Ω).(50)
As in Equation (42) the integrand in Equation (50) involves only (squared) scalar
products between the vectors !
vand !
wand the two vectors !b1and !b2of the basis. As
before this reflects the existence of the other two planes of reflection σwu and σvu as
indicated by the cartoon in Figure 1. We can again take recourse to the tensor identity
in Equation (43) and rewrite Equation (50) as
η=1
30dΩ[!b1·Q!
vv ·!b1+!b2·Q!
ww ·!b2
−!b2·Q!
vv ·!b2−!b1·Q!
ww ·!b1]α(Ω),(51)
where the tensors Q!
vv and Q!
ww are defined by complete analogy with the expression
given in Equation (36). We also emphasise that the expressionpresentedinEqua-
tion (51) agrees with Equation (19) in the paper of Camp and Allen [23].
That Equation (50) is a physically sensible definition of the biaxiality order parame-
ter can be established as follows. Consider the case in which αis isotropic with respect
to the axes !b1and !b2.Letthenxbe the cosine of the angles between each pair of the
four vectors !b1and !b2on the one hand and !
vand !
won the other hand. Thus, we are
dealing with four integrals of the form ±-1
−1dxx
2=±2
3which cancel so that η=0
in the absence of biaxial order.
Notice also, that even if αis isotropic with respect to !b1and !b2it may or may not
indicate alignment of the mesogens with respect to !b3so that S=0or S>0.Inthe
latter case, the C2-axis of the mesogens normal to the reflection plane σh(see Figure1)
would be strongly aligned while the mesogens can freely rotate about this axis. If they
can no longer rotate freely about the C2-axis we have to anticipate that S>0and
η>0.
We are now seeking to establish an alternative expression forηin terms of the Wigner
rotation matrices. To that end we remind the reader that because the mesogens are
centrosymmetric only even integers l,m,andncan arise. We also need to keep in
mind that we are deliberately neglecting any coupling between uniaxial and biaxial
degrees of freedom in the intermolecular interaction potential [see Equation (14)]. As
aconsequenceandaccordingtoMulder“...aspatiallyhomogeneousphasewillhavea
symmetry higher or equal to that of the constituent particles...” [21]. By that token,
α2
20 =α2
02 =0because these expansion coefficients refer to a coupling between uniaxial
and biaxial degrees of freedom [see Equations (A11) and (A12)]. We verified that in
our simulations this expectation is met.
Moreover, invoking the same symmetry considerations delineated in Section 2, it
seems obvious that we can alternatively consider
η=1
20dΩD2
±2±2(Ω)α(Ω)
=1
2$
lmn
αl
|m||n|s!
|m||n|0dΩD2
±2±2
∗(Ω)Dl
mn(Ω)(52)
17
where we employed Equation (31) and used D2
±2±2as shorthand notation for the ex-
pression D2
22 +D2
22 +D2
22+D2
22. Because of this definition D2
±2±2is real. Invoking now
the orthogonality of the Wigner rotation matrices [see Equation (21)] we obtain
η=1
2$
lmn
8π2
2l+1αl
|m||n|s!
|m||n|δl2δm±2δn±2(53)
which gives rise to four equivalent terms. Thus, using also Equation (32) together with
the definition of s!
mn (see above) we arrive at our final result
η=α2
22.(54)
Finally, using the definition of the Wigner rotation matricesgiveninEquation(5)
it is relatively straightforward to show that
D2
±2±2(Ω) = .cos2θ+1
/cos (2φ)cos(2χ)
−2cosθsin (2φ)sin(2χ).(55)
Thus, using this result in Equation (52) we obtain
η=0dΩ91
2.cos2θ+1
/cos (2φ)cos(2χ)
−cos θsin (2φ)sin(2χ):α(Ω)(56)
which agrees with the expression presented in Equation (17) in the paper by Camp
and Allen [23] and also matches Equation (3.9) in Mulder’s paper [21] as well as
Equation (5) in the paper by Straley [24].
4. Results
Throughout this work all physical quantities will be given inthecustomarydimen-
sionless (i.e., “reduced”) units. A comprehensive compilation of dimensionless units can
be found in Appendix B of the book by Allen and Tildesley [46]. The results pre-
sented below have been obtained in Monte Carlo (MC) simulations carried out in the
isothermal-isobaric ensemble. In this ensemble one generates a Markov chain with a lim-
iting distribution in configuration space proportional to exp{−β[U+PV −Nβ−1ln V]}
where Uis the total configurational potential energy, Pdenotes the pressure, and V
is the volume occupied by the N=5000mesogens.
Throughout this work, P=1.00.Oursimulationrunscomprise105MC steps con-
sisting of Nrandom rotations or displacements of the mesogens followed by one at-
tempt to change V.Thedecisionwhetheramesogen’scentreofmassistobedisplaced
or whether it is rotated about one of its three symmetry axes isdecidedatrandom
with equal on-average probability. Whether or not any of these attempts is accepted is
decided on the basis of a modified Metropolis criterion such that the aforementioned
limiting distribution in configuration space will eventually be reached. For details the
interested reader is referred to standard textbooks in the field [46, 47].
18
(a)
(b)
(c)
Figure 3. “Snapshots” of individual configurations obtained during the course of an MC simulation in the
isothermal isobaric ensemble. (a) Isotropic phase, T=1.120;(b)uniaxialnematicphase,T=1.030;(c)biaxial
nematic phase, T=0.894.Thesnapshotshavebeenobtainedforε!! =0.048.Seetextforanexplanationof
the colour code (see also Sections 2 and 4.1).
4.1. The formation of uniaxial and biaxial nematic phases
We begin the presentation of our results in Figure 3 by presenting “snapshots” of con-
figurations obtained from our MC simulations. We remind the reader that in reality we
are dealing with spherical particles (see discussion in Section 2) where the interaction
potential is only equipped with certain symmetry elements corresponding to molecules
that pertain to the point group D2h.Matchboxesarejustoneofseveralpossiblere-
alisations of objects belonging to that point group. To enhance the visualization of
our simulation-generated data we arbitrarily assigned the shape of matchboxes to the
mesogens under study.
The mesogens are coloured according to the following algorithm. Consider the three
largest eigenvalues λu
+,λv
+,andλw
+of the three tensors Quu,Qvv,andQww obtained
from Equation (37). These eigenvalues are sorted such that µ+=max(λu
+,λ
v
+,λ
w
+)
and µ−=min(λu
+,λ
v
+,λ
w
+);theremainingeigenvalueµ0is then taken to be equal to
that eigenvalue in the triplet λu
+,λv
+,andλw
+that is smaller than µ+but larger than
µ−.Theeigenvectorsassociatedwithµ+,µ0,andµ−are given by !
n+,!
n0,and!
n−,
respectively. They are also obtained as solutions of Equation (37).
We then introduce a triplet of real numbers bi=!
ui·!
n+,ri=!
vi·!
n−,andgi=;
wi·!
n0.
The triplet (ri,g
i,b
i)defines a point in the red, green, and blue (rgb) colour space.
Clearly, if the system is isotropic the mesogens appear in almost every colour of the
rgb colour space. This is the case in Figure 3(a). In a uniaxialnematicphase,biis
largest whereas riand giare nearly zero. As a result the mesogens are predominantly
coloured in shades of blue. This is the situation depicted in Figure 3(b). Finally, if
19
the system is in the biaxial nematic phase ri≃gi≃bi≃1.0which corresponds to a
light grey or white colour. From Figure 3(c) we therefore conclude that this snapshot
corresponds to a state of high biaxial order.
0.04 0.06 0.08 0.10
ε!!
0.90
1.00
1.10
1.20
1.30
1.40
T
0.00
0.20
0.40
0.60
0.80
1.00
Figure 4. Plot of the arithmetic mean order parameter O(defined in the text) [see Equations (42) and (50)] in
the plane spanned by the temperature Tand the dimensionless biaxial coupling constant ε!! [see Equation (14)].
The value of Ocan be read offthe attached colour bar. The white vertical lines correspond to temperature
scans for fixed biaxial coupling constants ε!! =0.03 [see Figure 5(a)], ε!! =0.05 [see Figure 5(b)], and ε!! =0.10
[see Figure 5(c).
To make this analysis more quantitative we consider in Figure4theorderparameter
O=1
2(S+η)corresponding to the arithmetic mean of the uniaxial and biaxial order
parameters. On account of its definition, O=0in the isotropic liquid phase. In a
perfectly ordered, uniaxial nematic phase O=0.5whereas 0.5≤O≤1.0if the
nematic phases possesses a certain degree of biaxial order. This is because in the
biaxial nematic phase both Sand ηare both fairly large.
The plot in Figure 4 shows that regardless of the value of ε!! only an isotropic liquid
phase is thermodynamically stable if the temperature Tis sufficiently high. As Tis
lowered the mesogens begin to form a nematic phase which is uniaxial at first but then
quickly becomes biaxial to a certain degree. The degree of biaxiality increases as Tis
lowered according to ones physical expectation.
As ε!! decreases the temperature of the onset of order decreases as well. This is
plausible because the formation of order is foiled by the kinetic energy the mesogens
possess on average. For sufficiently low values ε!! !0.04 it seems that even at the
lowest Tconsidered, biaxial order is completely absent or at least very small. The
overall topology of the phase diagram plotted in Figure 4 matches qualitatively the
one obtained by Sonnet et al. [5] at mean-field level (see their Figure 4).
These general features are illustrated in greater detail by the plots in Figure 5. For a
relatively small biaxial coupling constant ε!! =0.03 the plots in Figure 5(a) show that
as Tis lowered, both Sand ηare nearly zero as one would expect for an isotropic liquid.
At about T≃1.10,Sincreases rather abruptly and then keeps increasing slightly and
monotonically until T≃0.80 is reached. Because η≃0almost all the way down to
T≃0.80,theincreaseofSat about T≃1.10 signals the formation of a uniaxial
nematic phase. We cannot go to significantly lower temperatures T!0.80 because of
the onset of a formation of glassy phases.
For a somewhat larger biaxial coupling constant ε!! =0.05 we see from Figure 5(b)
that we have two characteristic temperatures instead of justone.Atbothtemperatures
20
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.90 1.00 1.10 1.20 1.30
S,η
(a)
S,η
(b)
S,η
T
(c)
Figure 5. Plots of the nematic order parameter S(•)andofthebiaxialorderparameterη(")asfunctions
of the temperature Tfor (a) ε!! =0.03,(b)ε!! =0.05,and(c)ε!! =0.10 (see Figure 4). The solid lines
are intended to guide the eye and the vertical dashed lines demarcate the estimated (a) isotropic-uniaxial
nematic (T≃1.10), (b) isotropic-uniaxial nematic (T≃1.14)anduniaxial-biaxialnematic(T≃0.99), and (c)
isotropic-biaxial nematic phase transitions (T≃1.24).
the order parameters increase markedly and independently. At the higher T≃1.14 ,
Sincreases notably whereas ηremains practically zero. Thus, at this temperature the
system undergoes a phase transition from an isotropic liquidtoauniaxialnematic
phase. As Tis lowered further, Skeeps increasing monotonically indicating that the
uniaxial nematic phase still becomes more ordered. At an evenlowerT≃1.00,the
biaxial order parameter ηstarts rising as well. This indicates that at this Twe have a
transition from a uniaxial to a biaxial nematic phase.
Finally, for the largest biaxial coupling constant ε!! =0.10 the plots in Figure 5(c)
reveal that we have only a single transition during which Sand ηsimultaneously rise
from about zero to fairly substantial values of S≃η≃0.75.Thus,atthecharacteristic
temperature T≃1.24 at which this happens we observe a transition directly from an
isotropic liquid to a biaxial nematic without an intermittent uniaxial nematic phase.
Notice also that the isotropic-biaxial nematic transition in Figure 5(c) occurs at the
highest temperature compared with Figures 5(a) and 5(b). This is in line with one’s
physical intuition because the formation of a biaxial nematic is controlled by the largest
biaxial coupling constant considered in the three parts of Figure 5. Therefore ordered
phases may form even though the thermal energy of the mesogensislargest.
21
(a)
(b)
b1b2
b3
b1b2
b3
Figure 6. Sketch of the equatorial planes introduced in Section 3.2 to promote the analysis of the orientation
distribution functions αuand αv.ThecolourcodeisthesameastheoneintroducedinFigure2.The rings
represent the respective equators in which αuand αvare computed after the rotations of the eigensystem
described in Section 3.2. Arrows in the equatorial planes indicate preferred orientations of the mesogens if the
system exhibits perfect uniaxial (a) or biaxial order (b).
4.2. Orientational and positional structure
−
π
20π
2
˜ϕu
0
π
2
π
θu
(a)
−
π
20π
2
˜ϕu
0.0
0.7
1.4
(b)
Figure 7. Plots of the orientation distribution function αuin the "ϕu–θuplane where "ϕu=ϕu−π(ϕu∈
[0,2π])hasbeenshiftedtocentertheplotson"ϕ=0.Thevalueofαucan be read offthe attached colour bar.
In all cases ε!! =0.048;(a)T=1.04,(b)T=0.91.
After focusing on the formation of uniaxial and biaxial nematic phases in the pre-
ceding section we now elucidate details of the structure of these phases. Before showing
any results we find it useful to remind the reader that we compute the odf’s not with
the original set of vectors {!
ui}and {!
vi}but with the rotated vectors {!
u!
i}and {!
v!
i}in-
stead. The rotated vectors are obtained from Equations (38) as explained in Section 4.
Hence, vectors of the sets {!
u!
i}and {!
v!
i}lie in their respective equatorial planes.
In a perfectly ordered uniaxial nematic phase all vectors of the set {!
u!
i}are pointing
in the direction of the axis ±!b2of the equatorial plane represented by the red ring in
Figure 6(a). Because the perfect nematic phase is uniaxial thevectorsoftheset{!
v!
i}
are distributed uniformly along the equatorial plane represented by the blue ring in
Figure 6(a). Consequently, in a uniaxial nematic phase that is not perfectly ordered one
would expect αuto be centered on the axis labelled !b2in Figure 6(a) and represented
by the red arrow. The corresponding odf αvwould be uniform along the equator in the
22
plane !b1–!b2with a certain width in the ±!b3-direction.
In the perfectly ordered biaxial nematic phase, schematically represented by the plot
in Figure 6(b), the situation remains unaltered as far as αuis concerned. However, in
aperfectbiaxialnematicphasethevectorsoftheset{!
v!
i}are all pointing along the
±!b1-axis which is chosen at random with an equal probability. Forabiaxialphase
that is not perfectly ordered one expects an odf consisting offourmoreorlesssharp
maxima in the !b1–!b2plane separated by an angle of π.
−
π
20π
2
˜ϕv
0
π
2
π
θv
(a)
−
π
20π
2
˜ϕv
0.0
0.4
0.8
(b)
Figure 8. As Figure 7 but for αv.
We begin our discussion by presenting in Figure 7 plots of αu.Weseethatinboth
parts of the figure a single centrosymmetric spot appears thatturnsouttobehigher
at the lower T.Thespotsarecenteredonθu=π
2and <ϕu=0.Tomakecontactwith
the cartoon shown in Figure 6(a) we see that the axis !b1should be orthogonal to the
paper plane. We emphasise that αuexhibits a second maximum at <ϕu=πwhich we
do not show in Figure 7. Because the spots in Figure 7 are centrosymmetric we can,
however, exclude the existence of any extrinsic biaxiality.
We can also interpret the plot in Figure 7(a) in a slightly different way. It is clear
that αureflects the distribution of the orientations {!
ui}about the nematic director !
n.
Therefore, !
nis pointing in a direction normal to the paper plane. The centrosymmetric
spot of the odf in Figure 7(a) indicates that !
nis a C∞rotation axis.
Unfortunately, we cannot distinguish uniaxial from biaxialnematicphasesonthe
basis of αualone which holds for the presence of intrinsic biaxiality. In order to dis-
criminate between both symmetries we need to consider the odfcalculatedforasecond
symmetry axis. To that end we consider in Figure 8 plots of αv.Thedataplottedin
Figure 8(a) can be represented by a band-like structure of a certain width around
θv=π
2extending along the <ϕv-axis. The existence of the band reflects the uniform
distribution of {!
v!
i}along the equatorial plane shown in Figure 6(a). The width of
the band is a measure of the average deviation of the vectors {!
v!
i}from the !b3-axis.
Thus, on account of the combined information contained in Figures 7(a) and 8(a) we
conclude that the system is in the uniaxial nematic phase.
At the lower T=0.91 the corresponding plot presented in Figure 8(b) now exhibits
asequenceofisolated,slightlyelongatedspotswherethemajor axis is parallel to the
23
0.0
1.0
2.0
3.0
4.0
1.0 2.0 3.0 4.0
g(r12)
r12
Figure 9. Plot of the radial pair correlation function g(r12)as a function of the centre-of-mass distance r12
for a pair of mesogens at a temperature T=0.89 and ε!! =0.05 at which the system is in the biaxial nematic
phase (see Figure 4).
<ϕv-axis and the minor axis is parallel to the θv-axis. Because of the arguments put
forth when discussing the cartoon shown in Figure 6(b), thesecentresofthespotsare
separated by an angle increment of ∆<ϕv=π
2.Thus,weconcludethattheplotsshown
in Figures 7(b) and 8(b) are indicative of a biaxial nematic phase. Notice again that
we are showing only one of the two maxima of which αvconsists in the biaxial nematic
phase.
In other words, because the centrosymmetric character of thespotofαuin Fig-
ure 7(b) remains unchanged compared with Figure 7(a), !
nremains a C∞rotation axis
even if the system is biaxial. By the same token, αvplotted in Figure 8(b) reflects the
distribution of orientations of the mesogens about the secondary symmetry axis which
is also pointing in a direction normal to the paper plane. However, the difference is
that this secondary symmetry axis has C2character unlike !
n.Weascribethistothe
lack of coupling between uniaxial and biaxial degrees of freedom in our model.
We now turn to an analysis of orientational and positional correlations in our system.
To make sure that our system is still fluidic we consider the radial pair correlation
function defined as [38]
g(r12)= 1
(8π2)200dΩ1dΩ2g(r12,Ω1,Ω2).(57)
AplotofthisquantityispresentedinFigure9forathermodynamic state point at
which the system is in a stable biaxial nematic phase. As one can see in Figure 9,
spatial correlations are still short-range and decay withinacoupleofdiametersσof the
spherical core of our matchbox-symmetric mesogens. However, there is a noticeable shift
of the maxima to progressively smaller values of r12 that is untypical for a simple fluid.
Whereas the first peak of gappears at about the value at which uLJ in Equation (13)
passes through its minimum, higher-order peaks appear at non-integer multiples of the
position of the first maximum. This feature is ascribed to the different (non-spherical)
symmetry of the biaxial nematic phase. However, we emphasizethatthecurveshown
in Figure 9 is free of any pathological features that would indicate an artificial (and
unwanted) glassy structure.
We are now turning to the orientation correlation functions.Theseareaccessible
through Equations (33). Generally speaking, the orientation correlation functions plot-
24
0.0
0.2
0.4
S2
,2η2
0.0
0.4
0.8
2η2
S2
0.5
1.0
1.5
1.0 2.0 3.0 4.0
S2
2η2
!Sn1n2
220 "shell
(a)
!Sn1n2
220 "shell
(b)
!Sn1n2
220 "shell
r12
(c)
Figure 10. Plots of the orientation correlation functions #Sn1n2
220 $shell as functions of the distance r12 between
apairofmesogens;(
")n1=n2=0,(•)n1=±2,n2=±2.Inallpartsofthefigure,ε!! =0.05;(a)T=1.12,
(b) T=1.03,(c)T=0.89.AlsoindicatedarethelimitingvaluesofS2and 2η2obtained from Equations (42)
and (50). Note that in all three parts of the figure, the limiting values S2and 2η2are approached from above
as r12 increases according to one’s physical intuition.
ted in the three parts of Figure 10 exhibit far less structure than the radial pair cor-
relation function shown in Figure 9. In the isotropic phase, the plots presented in Fig-
ure 10(a) indicate that )Sn1n2
220 *shell decay monotonically and vanish completely within
acoupleofdiametersσof the spherical cores. It turns out that 1S±2±2
220 2shell vanishes
more rapidly than 1S00
2202shell indicating that biaxial short-range order is lost earlier.
The data shown in Figure 10(b) pertain to a state point at whichthesystemisina
uniaxial nematic phase. This is concluded because 1S±2±2
220 2shell is still short-range and
vanishes on a molecular length scale whereas 1S00
2202shell has become long-range with
no appreciable tendency to decay towards zero.
Finally, in the biaxial nematic phase 1S00
2202shell as well as 1S±2±2
220 2shell exhibit long-
range order as one might have expected. However, it is particularly gratifying that in
all three parts of Figure 10 the curves decay to the values predicted by Equations (35a)
and (35b). This lends credibility to our simulations becausethelimitingvaluesofS2
and 2η2were computed independently from Sand ηas given by Equations (42) and
(50).
In view of the results just discussed it seems instructive to focus on the decay of
25
-8.0
-7.0
-6.0
-5.0
-4.0
-3.0
-2.0
1.0 1.2 1.4 1.6 1.8
ln ∆Sn1n2
220
ln r12
Figure 11. Double-logarithmic plots of ∆Sn1n2
220 as functions of r12 ;(•)n1=n2=0,(")n1=±2,
n2=±2.Thestraightfinesarelinearfitstothediscretedatapoints.ThedatashowncorrespondtoT=1.08
and ε!! =0.05
short-range orientational correlations. To that end we consider
∆S00
220 ≡1S00
2202shell −S2(58a)
∆S±2±2
220 ≡1S±2±2
220 2shell −2η2(58b)
which vanish by definition in the limit of r12 →∞.Moreover,letusassumethat
both quantities defined in Equations (58) decay exponentially. Thus, we can define a
correlation length ξn1n2as the slope of the straight lines shown in Figure 11.
Repeating the analysis illustrated by the plots in Figure 11 for different temperatures
we obtain the curves shown in Figure 12. The plots in Figure 12(a) show that ξ00
decays towards the temperature at which a uniaxial nematic phase forms; it assumes a
constant value of about 3.0for all lower T.Theplotofξ±2±2increases monotonically
throughout the temperature range considered. This correlates nicely with the absence
of a uniaxial-biaxial nematic transition but a steady increase of short-range biaxial
correlations in the corresponding Figure 5(a). It seems plausible that ξ00 assumes a
constant value once the uniaxial nematic phase has formed because of our definition
of the correlation lengths.
Similarly, the plots in Figure 12(b) show very similar trends. Towards the isotropic-
uniaxial nematic transition, ξ00 decays and then assumes a constant value for all lower
T.Onthecontrary,ξ±2±2increases towards the uniaxial-biaxial nematic transition
and then decays again for all lower T.Notice,inparticular,thatbelowtheuniaxial-
biaxial-nematic transition ξ00 ≈ξ±2±2.Thelattereffectisseenmoreclearlyintheplots
presented in Figure 12(c) where no noticeable variation of ξ00 and ξ±2±2is detected
by lowering Tbelow the isotropic-biaxial nematic transition temperature at T≃1.24.
5. Discussion and conclusions
By means of MC simulations in the isothermal-isobaric ensemble we investigate the
formation and properties of uniaxial and biaxial nematic phases in a thermotropic
liquid crystal. The liquid crystal is composed of matchbox-symmetric mesogens which
26
0.00
2.00
4.00
6.00
0.00
2.00
4.00
6.00
0.00
2.00
4.00
6.00
0.90 1.00 1.10 1.20 1.30
ξn1n2
(a)
ξn1n2
(b)
ξn1n2
T
(c)
Figure 12. Plots of the correlation length ξn1n2(see text) as a function of the temperature Tfor n1=n2=0
(•)andn1=±2,n2=±2("). Solid lines are fits to the discrete data points intended to guide the eye; (a)
ε!! =0.03,(b)ε!! =0.05,(c)ε!! =0.10.Verticaldashedlinesindicatethetemperaturesatwhichordered
phases form (cf., Figure 5).
possess an intrinsic biaxiality. This is reflected by the factthatthemesogenshave
amajorandtwoequivalentminorsymmetryaxesandthuspossess three mutually
orthogonal reflection planes.
In a number of recent computer simulation studies [28, 48–50]theformationof
ordered phases was studied for models in which the molecular constituents interact
via purely hard repulsive potentials. In such systems it turns out to be a bit difficult
to form and stabilise ordered structures. This could be because of the purely entropic
nature of the interactions between the particles which may betooweak.Inaddition,
any attempt to change the orientation of the molecules duringthecourseofasimulation
is difficult to accomplish. One way to overcome this difficulty istouseexternalfields
[50]. Another disadvantage of such model systems in computersimulationsisthatthe
interaction strength cannot be controlled easily.
We are proposing a different strategy in this work. Rather thanmodelingthe
shape of the mesogens explicitly we take spherical molecules and endow them with an
orientation-dependent attraction that accounts for the symmetry elements of molecules
of the desired point group. This philosophy follows in spirittheearlierworkofMaier
and Saupe [41]. Thus, in a sense we are dealing with pseudo-shaped mesogens.
Systems composed of pseudo-shaped mesogens generally allowforanequilibration
that is orders of magnitude faster than for models with explicit shape dependence.
At the same time pseudo-shaped mesogens do exhibit physically sensible properties
27
as we demonstrated in a number of earlier publications [51–53]. From our point of
view purely entropic interactions have quite a few shortcomings as far as simulation
studies are concerned. However, we stress that their use can be quite advantageous in
theoretical approaches such as classical density functional theory [54] on account of
the vanishing range of the entropic interactions.
In this work the interaction between a pair of mesogens is described by a model
potential where the orientation dependence is the same as that proposed originally
by Luckhurst and Romano [30]. The model proposed by these authors has also been
used later by Longa and Pająk [31]. In the Luckhurst-Romano model the orientation
dependence of the interactions follows from an expansion of the potential in terms of
Stone’s rotational invariants %Sn1n2
l1l2l&[29]. These functions are rotationally invariant
products of the Wigner rotation matrices Dl
mn [38, 39]. However, we are using an
off-lattice model unlike the Luckhurst-Romano model.
We emphasise that the expansion in Wigner rotation matrices of the interaction
potential and of the orientation-dependent pair correlation function can also be cast
in terms of the contraction of Cartesian tensors. These tensors describe the coupling
between uniaxial degrees of freedom on the one hand and biaxial degrees of freedom on
the other hand. The latter formulation is identical with the one proposed by Sonnet et
al. [5] and by Sonnet and Virga [36]. We have verified here that bothformulationsare
indeed the same. The proof of this equivalence rests upon a ladder-operator algebra
[34] borrowed from the quantum-mechanical theory of angularmomentum[35].
In summary, we are presenting three key results in this work. The first of these
concerns the good agreement of the topology of our phase diagram (see Figure 4) with
the one presented by Sonnet et al. [5] and obtained within their mean-field treatment.
Within the parameter space investigated here three types of phase transitions can
occur. The first of these exhibits only an isotropic to uniaxial nematic transition. The
second one allows for a sequence of separate isotropic-uniaxial nematic followed by
auniaxial-biaxialnematictransition.Thethirdscenariotakes one directly from an
isotropic to a biaxial nematic phase. Relatively recent experimental evidence for the
latter two scenarios seems to be provided by 13CNMRspectroscopy[13].Theoretically,
the above sequence of phase transitions is predicted by the phenomenological Ginzburg-
Landau theory [31].
The second key result concerns an analysis of the odf. From a separate analysis
of the odf’s referring to the respective major and (one of the two) minor symmetry
axes we can unambiguously distinguish between uniaxial nematic and biaxial nematic
orientational structures. At its core, the analysis is enabled by properly rotating the
molecular major and minor symmetry axes such that they alwayslieintheequatorial
plane of the unit sphere. This way we are able to combine an enhanced resolution of
the odf with a maximum statistical accuracy with which the odf’s are accessible in a
computer simulation. The idea to employ a rotating rather than a fixed reference frame
has been proposed in the work of Skutnik et al. [7].
On the basis of the odf one can unambiguously distinguish between a system with
intrinsic and one with extrinsic biaxiality. For the latter it was recently demonstrated
[7] that in a binary mixture of two uniaxial components a globally biaxial nematic phase
can form if the orientation dependence of the interaction between two unlike mesogens
is tuned such that a T-shaped arrangement is energetically favoured. In the biaxial
phase the odf is elliptically deformed in the <ϕ–θplane where the two major axes of the
ellipses are orthogonal to one another. In a system with extrinsic biaxiality only one
odf conveys the entire information about the structure of thebiaxialnematicphase.
28
Because of the reflection plane σhthe αuexhibits two equivalent circular maxima that
are separated by an angle increment ∆<ϕu=πregardless of whether one is in the
uniaxial or biaxial nematic phase.
In the uniaxial nematic phase αvconsists of a strip centered on θu=π
2because there
is no particular order with respect to the two minor symmetry axes. In the biaxial
nematic phase, the reflection plane σwu
vsuggests that the odf αvshould also exhibit
two maxima separated by an angle ∆<ϕv=π.However,oneobservestwomaximathat
are but separated by ∆<ϕv=π.Thisisbecausethetwominormolecularaxes!
viand !
wi
are degenerate because we neglect the coupling between uniaxial and biaxial degrees
of freedom in our model. Thus, our mesogens pertain to the D2h point group.
The simple structure of the odf’s reflects the overall symmetry of the interaction
potential that does not contain a coupling between uniaxial and biaxial degrees of
freedom. Therefore, the two maxima of αuas well as the two maxima of αvin the
biaxial nematic phase have equal height. It is anticipated that when this symmetry is
broken by allowing for this coupling, the centrosymmetric spots of αuwill be elliptically
deformed as well. In other words, the C∞of !
nwill presumably change to the lower C2
symmetry already seen in αvin the present work.
The difference between that earlier work and the present one isthathereweare
employing as a reference system the system of eigenvectors ofthethreealignment
tensors referring to the major and the two minor symmetry axesofthemesogens.
This approach follows in spirit the one used by Allen [22] and by Camp and Allen [23]
in their calculation of the uniaxial and biaxial order parameters. These authors did,
however, make no attempt to analyse the odf.
The third key result is related to our analysis of the orientation-dependent pair cor-
relation function when their expansion is made in terms of Stone’s rotational invariants
[29]. This expansion was proposed earlier by Blum and Torruella [55] who did, however,
not show any actual data for the expansion coefficients. We demonstrate in this work
that the expansion coefficients become long-range when an ordered phase has formed.
As we show here the limiting values of the correlation functions in general involve a
subtle coupling between different order parameters. The simpler result obtained quite
some time ago by Zannoni [37, 43] is recovered as a special caseofourtreatment.
Based upon the assumption that the decay of the expansion coefficients towards their
long-range value is exponential, we can extract a correlation length for the short-range
correlations. This correlation length assumes a constant value once the “macroscopi-
cally” ordered phase has formed. For thermodynamic states prior to the formation of
“macroscopically” ordered phases the temperature dependence of the correlation length
signals uniaxial or biaxial “pre-ordering”. The biaxial pre-ordering is a precursor of a
biaxial nematic phase eventually triggering the latter. This situation is very much akin
to that encountered experimentally where first tetrapodes are being synthesised that
are of biaxial symmetry only on a local length scale. These then initiate the formation
of a globally biaxial nematic phase given the right thermodynamic conditions [17].
Last but not least, we emphasise that the approach employed here is flexible enough
to be extended to treat interaction potentials in models of bent-core mesogens. These
mesogens have two rather than three reflection planes and thusasomewhatlower
symmetry compared with the matchbox-symmetric mesogens studied here. In these
systems undulated biaxial orientational order can emerge [see Figure 1(e) of Ref. 56].
Such undulated structures are conceivable if for the presentmodeltheuniaxial-and
biaxial-tensor contributions are modified to allow for a certain handedness of the inter-
action potential in two directions. That way it seems conceivable that one could study
systems exhibiting chiral C∗phases. Work along these lines is currently in progress.
29
Appendix A. Cartesian tensors
The purpose of this Appendix is to derive expressions for Sn1n2
220 arising in the expansion
of uaniso in Equation (6). We begin the discussion with S00
220 expressing it as a Cartesian
tensor according to [29]
2√5S00
220 (Ω1,Ω2,Ω12)=3(
!
u1·!
u2)2−1,(A1)
where !
uiis given in Equation (11a). Using the identity (see Section 2)
!
u1!
u1:!
u2!
u2=Tr(
!
u1!
u1·!
u2!
u2)=(
!
u1·!
u2)2(A2)
and the definition of Ugiven in Equation (10a), it is a simple matter to verify that
S00
220 (Ω1,Ω2,Ω12)= 1
√5
3
2U1:U2(A3)
To proceed it turns out to be most convenient [29] to draw upon the angular-
momentum algebra in quantum mechanics [35] and introduce thedifferentialoperator
Lj=−i∇Ωj(A4)
such that [29]
Lj(!
vj!
wj!
uj)=
0−i!
uji!
wj
i!
uj0−i!
vj
−i!
wji!
vj0
.(A5)
Again by complete analogy with quantum mechanics we introduce the ladder operators
Lj±=Ljv ±iLjw (A6)
satisfying the equation
L1±Sn1n2
l1l2l(Ω1,Ω2,Ω12)
=[l1(l1+1)−n1(n1∓1)]1/2S(n1∓1)n2
l1l2l(Ω1,Ω2,Ω12).(A7)
The subscript jrefers to the mesogen the ladder operator is acting upon.
Applying Equation (A7) to S00
220 and using Table 2 of the paper by Stone [29] we
obtain
L1+S00
220 =√6S10
220
=1
2√5(L1v+iL1w)33(!
u1·!
u2)2−14
=3
√5(!
u1·!
u2)(!
v1+i!
w1)·!
u2(A8)
where we dropped the arguments of S00
220 and S10
220 for the time being and also used
30
Equation (A5). From the previous expression it then follows that
S10
220 =53
10 (!
u1·!
u2)(!
v1+i!
w1)·!
u2.(A9)
Applying L1+ one more time to this result we obtain
L1+S10
220 =2S20
220
=53
10L1+ (!
u1·!
u2)(!
v1+i!
w1)·!
u2
=53
10 [(!
v1+i!
w1)·!
u2]2
=53
10 3(!
v1·!
u2)2−(!
w1·!
u2)2+2i(!
v1·!
u2)(!
w1·!
u2)4.(A10)
We can finally solve this expression for S20
220 and obtain
S20
220 =1
253
10 [B1:U2+2i(!
v1·!
u2)(!
w1·!
u2)] ,(A11)
where we used Equations (10a), (10b), (A2), and the fact that the tensors !
u2!
u2,!
v1!
v1,
and !
w1!
w1satisfy the distributive law. Notice also that we do not have to worry whether
or not we have to take the transpose of any tensor arising in theabovemanipulations
because they are all real and symmetric.
The derivation of S20
220 proceeds exactly the same way. The calculation can be
shortened considerably by noting that in successive applications of L1−,L1−!
u1=
−(!
v1−i!
w1)whereas L1−(!
v1−i!
w1)=0.Thus,attheendofthedayweobtain
S20
220 =1
253
10 [B1:U2−2i(!
v1·!
u2)(!
w1·!
u2)] (A12)
which is identical with the expression in Equation (A11) except for the sign of the
imaginary term. Thus, because of the symmetry considerations delineated in Section 2
it seems sensible to introduce
S20
220 +S20
220 =53
10B1:U2.(A13)
Finally, repeating the above analysis with L2±instead of L1±we would finally arrive
at
S02
220 +S02
220 =53
10U1:B2(A14)
and thus the coupling between uniaxial and biaxial degrees offreedomtouaniso would
be proportional to B1:U2+U1:B2.However,intheinterestofminimumcomplexity
of the model we shall follow Sonnet et al. [5] and disregard this type of coupling
altogether. This last expression can be obtained by symmetry arguments. Because the
31
mesogens are identical we can interchange the labels 1and 2and replace n1by n2and
vice versa.
Finally, using the same calculus as above and successively applying L2±twice to
Equations (A11) and (A12) we obtain four equations. For symmetry reasons discussed
in Section 2 we can use Equation (9). One notices that each termontheright-hand
side of this latter expression is of the general form
1
451
5[(!
v1∓i!
w1)·(!
v2∓i!
w2)]2,(A15)
where the plus sign applies if the corresponding integer of Sn1n2
220 on the right-hand side
of Equation (9) is negative and vice versa. Using the definition of the biaxial tensor B
given in Equation (10b) as well as Equation (B.16) of Ref. 38 we obtain
S±2±2
220 =51
5B1:B2(A16)
as our final result.
Appendix B. Limiting value of the orientation correlation function
In this Appendix we rationalise Equations (33) for the limiting (r12 →∞)behaviour
of the orientation correlation functions introduced in Equations (26a) and (26b). Using
Equations (30) and (31) we have
1Sn1n2
ll02∞=(−)l
√2l+1"2l+1
8π2#2$
m
(−)m
×$
l1m1k1
s|m1||k1|αl1
|m1||k1|
×0dΩ1Dl1
m1k1(Ω1)Dl
mn1
∗(Ω2)
×$
l2m2k2
s|m2||k2|αl2
|m2||k2|
×0dΩ2Dl2
m2k2(Ω2)Dl
mn2
∗(Ω2),(B1)
where the extra factors in the summations over l1,2,m1,2and k1,2are introduced for
the same reason as in Equation (18). Using Equation (A.93) of Ref. 38 this can be
rewritten as
1Sn1n2
ll02∞=(−)l
√2l+1$
m
(−)m$
l1m1k1$
l2m2k2
s|m1||k1|s|m2||k2|
×αl1
|m1||k1|αl2
|m2||k2|δll1δll2δn1k1δn2k2δm1mδm2m.
(B2)
32
We restrict the evaluation of this expression to even lfor reasons already stated in
Section 2. In addition, the Kronecker symbols cause several of the summations to
collapse so that we can rewrite the previous equation more compactly as
√2l+11Sn1n2
ll02∞=$
m1m2m
(−)ms|m1||n1|s|m2||n2|
×αl
|m1||n1|αl
|m2||n2|δm1mδm2m.(B3)
Noticing that the integers m1and m2have to be even on account of the symmetry of
the odf [21] and evaluating the last two Kronecker symbols in the previous expression
we obtain
√2l+11Sn1n2
ll02∞=
l
$
m=0
αl
m|n1|αl
m|n2|2δn10+δn20−2,(B4)
where we evaluated s|m1||k1|s|m2||k2|using Equation (32). Notice that in Equation (B4),
the summation is carried out only for positive mbecause terms for positive and negative
values of mare identical. This brings in an extra factor of 2except for m=0.Tocorrect
for this we have introduced 2−δm0.Inaddition,itshouldbeborneinmindthaton
account of the symmetry of the odf, αl
|m1||k1|and αl
|m2||k2|vanish for odd values of m1,
m2,k1,ork2from the outset. From Equation (B4) it is easy to derive Equations (33)
and (34).
References
[1] P.G. de Gennes and J. Prost, The physics of liquid crystals,2nded.(OxfordUniversity
Press, Oxford, 1995).
[2] M.J. Freiser, Phys. Rev. Lett. 24, 1041–1043 (1970).
[3] R. Alben, J. Chem. Phys. 59, 4299–4304 (1973).
[4] J.H. Lee, T.K. Lim, W.T. Kim and J.I. Jin, J. Appl. Phy. 101, 034105 (2007).
[5] A.M. Sonnet, E.G. Virga and G.E. Durand, Phys. Rev. E 67, 061701 (2003).
[6] J.C. Eichler, R.A. Skutnik, A. Sengupta, M.G. Mazza and M.Schoen,Mol.Phys.117,
3715–3733 (2019).
[7] R.A. Skutnik, L. Lehmann, S. Püschel-Schlotthauer, G. Jackson and M. Schoen, Mol.
Phys. 117, 2830–2845 (2019).
[8] A. Cuetos, A. Galindo and G. Jackson, Phys. Rev. Lett. 101, 237802 (2008).
[9] L.J. Yu and A. Saupe, Phys. Rev. Lett. 45, 1000–1003 (1980).
[10] E.A. Oliveira, L. Liebert and A.M. Figueiredo Neto, Liquid Crystals 5, 1669–1675 (1989).
[11] A. Stroobants and H.N.W. Lekkerkerker, J. Phys. Chem. 88, 3669–3674 (1984).
[12] H.H. Wensink, G.J. Vroege and H.N.W. Lekkerkerker, Phys. Rev. E 66, 041704 (2002).
[13] R.Y. Dong, S. Kumar, V. Prasad and J. Zhang, Chem. Phys. Lett. 448, 54–60 (2007).
[14] M. Lehmann, C. Köhn, J.L. Figueirinhas, G. Feio, C. Cruz and R.Y. Dong, Chem. Eur.
J. 16, 8275–8279 (2010).
[15] B.R. Acharya, A. Primak and S. Kumar, Liquid Crystals Today 13, 1–4 (2004).
[16] E. van den Pol, A.V. Petukhov, D.M.E. Thies-Weesie, D.V.ByelovandG.J.Vroege,Phys.
Rev. Lett. 103, 258301 (2009).
[17] C. Cruz, J.L. Figueirinhas, D. Filip, G. Feio, A.C. Ribeiro, Y. Frère, T. Meyer and G.H.
Mehl, Phys. Rev. E 78, 051702 (2008).
[18] A. Eremin, S. Diele, G. Pelzl, H. Nádasi, W. Weissflog, J. Salfetnikova and H. Kresse,
Phys. Rev. E 64, 051707 (2001).
33
[19] D.Y. Kim, L. Wang, Y. Cao, X. Yu, S.Z.D. Cheng, S.W. Kuo, D.H. Song, S.H. Lee, M.H.
Lee and K.U. Jeong, J. Mater. Chem. 22, 16382–16389 (2012).
[20] B. Tjipto-Margo and G.T. Evans, J. Chem. Phys. 94, 4546–4556 (1991).
[21] B.M. Mulder, Liq. Cryst. 1, 539–551 (1986).
[22] M.P. Allen, Liq. Cryst. 8, 499–511 (1990).
[23] P.J. Camp and M.P. Allen, J. Chem. Phys. 106, 6681–6688 (1998).
[24] J.P. Straley, Phys. Rev. A 10, 1881 (1974).
[25] R. Berardi, C. Fava and C. Zannoni, Chem. Phys. Lett. 236, 462–468 (1995).
[26] R. Berardi and C. Zannoni, J. Chem. Phys. 113, 5971–5979 (2000).
[27] F. Biscarini, C. Chiccoli, P. Pasini, F. Semeria and C. Zannoni, Phys. Rev. Lett. 75, 1803
(1995).
[28] A. Cuetos, M. Dennison, A. Masters and A. Patti, Soft Matter 13, 4720–4732 (2017).
[29] A.J. Stone, in The molecular physics of liquid crystals,editedbyG.R.Luckhurstand
G. W. Gray (Academic Press, London, 1979), Chap. 2, pp. 31–50.
[30] G.R. Luckhurst and S. Romano, Mol. Phys. 40, 129–139 (1980).
[31] L. Longa and G. Pająk, Liq. Cryst. 32, 1409–1417 (2005).
[32] L. Longa, P. Grzybowski, S. Romano and E. Virga, Phys. Rev. E 71, 051714 (2005).
[33] Y. Martínez-Ratón, S. Varga and E. Velasco, Phys. Chem. Chem. Phys. 13, 13247–13254
(2011).
[34] A.J. Stone, Mol. Phys. 36, 241–256 (1978).
[35] A. Messiah, Quantum mechanics,Vol.2(NorthHolland,Amsterdam,1969).
[36] A.M. Sonnet and E.G. Virga, Dissipative ordered fluids: theories for liquid crystals
(Springer Science & Business Media, New York, 2012).
[37] C. Zannoni, in The molecular physics of liquid crystals,editedbyG.R.Luckhurstand
G. W. Gray (Academic Press, London, 1979), Chap. 3, pp. 51–83.
[38] C.G. Gray and K.E. Gubbins, Theory of molecular fluids,Vol.1(ClarendonPress,Oxford,
1984).
[39] D.M. Brink and G.R. Satchler, Angular momentum,2nded.(ClarendonPress,Oxford,
1968), pp. 146–148.
[40] A.A. Albert and B. Muckenhoupt, Mich. Math. J. 4, 1–3 (1957).
[41] W. Maier and A. Saupe, Z. Naturforsch. A 15, 287–292 (1960).
[42] M. Schoen, A.J. Haslam and G. Jackson, Langmuir 33, 11345–11365 (2017).
[43] C. Zannoni, in The molecular physics of liquid crystals,editedbyG.R.Luckhurstand
G. W. Gray (Academic Press, London, 1979), Chap. 9, pp. 191–220.
[44] W.B. Streett and D. Tildesley, Proc. Roy. Soc. (London) A348, 485–510 (1976).
[45] R. Eppenga and D. Frenkel, Mol. Phys. 52, 1303–1334 (1984).
[46] M.P. Allen and D.J. Tildesley, Computer simulation of liquids,2nded.(Oxforduniversity
press, Oxford, 2017).
[47] D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to appli-
cations,Vol.1,2nded.(AcademicPress,SanDiego,2001).
[48] B.S. John and F.A. Escobedo, J. Phys. Chem. B 109, 23008–23015 (2005).
[49] B.S. John, C. Juhlin and F.A. Escobedo, J. Chem. Phys. 128, 044909 (2008).
[50] A. Cuetos, E.M. Rafael, D. Corbett and A. Patti, Soft Matter 15, 1922–1926 (2019).
[51] S. Schlotthauer, R.A. Skutnik, T. Stieger and M. Schoen,J.Chem.Phys.142, 194704
(2015).
[52] S. Püschel-Schlotthauer, T. Stieger, M. Melle, M.G. Mazza and M. Schoen, Soft Matter
12, 469–480 (2016).
[53] S. Püschel-Schlotthauer, V.M. Turrión, T. T. Stieger, R. Grotjahn, C. Hall, M.G. Mazza
and M. Schoen, J. Chem. Phys. 145, 164903 (2016).
[54] M. Schoen and C. Hoheisel, Mol. Phys. 57, 65–79 (1986).
[55] L. Blum and A.J. Torruella, J. Chem. Phys. 56, 303–310 (1972).
[56] L.A. Madsen, T.J. Dingemans, M. Nakata and E.T. Samulski, Phys. Rev. Lett. 92, 145505
(2004).
34
