Surprisingly long-ranged intermolecular correlations begin to appear in isotropic (orientationally disordered) phases of liquid crystal forming molecules when the temperature or density starts to close in on the boundary with the nematic (ordered) phase. Indeed, the presence of slowly relaxing, strongly orientationally correlated, sets of molecules under putatively disordered conditions (“pseudo-nematic domains”) has been apparent for some time from light-scattering and optical-Kerr experiments. Still, a fully microscopic characterization of these domains has been lacking. We illustrate in this paper how pseudo-nematic domains can be studied in even relatively small computer simulations by looking for order-parameter tensor fluctuations much larger than one would expect from random matrix theory. To develop this idea, we show that random matrix theory offers an exact description of how the probability distribution for liquid-crystal order parameter tensors converges to its macroscopic-system limit. We then illustrate how domain properties can be inferred from finite-size-induced deviations from these random matrix predictions. A straightforward generalization of time-independent random matrix theory also allows us to prove that the analogous random matrix predictions for the time dependence of the order-parameter tensor are similarly exact in the macroscopic limit, and that relaxation behavior of the domains can be seen in the breakdown of the finite-size scaling required by that random-matrix theory.

Liquid crystals are the quintessential example of liquids whose order spans a significant range of possibilities. In the simplest examples, liquids switch between an isotropic distribution of molecular orientations and a uni-axial nematic distribution (one in which there is a single preferred axis for molecular alignment) as the temperature, density, or concentrations are varied.1 What we are mostly concerned with in this paper is what we can say about the nature of whatever orientational order there is in the isotropic phase of such systems.

The designation “isotropic” means that there is no macroscopically long-ranged order, and the first-order character of the isotropic/nematic transition2 would seem to suggest that an isotropic phase should remain placidly unaware of the possibility of nematic order until the transition is reached. Yet, as has been shown by light scattering,3–7 video microscopy,8,9 and a variety of ultrafast spectroscopic methods,10–17 orientational correlation lengths and times both grow rapidly as the transition is approached from the isotropic side.2,18–23 Landau theories of the isotropic/nematic transition,2,24 in fact, predict that there ought to be a hidden isotropic/nematic critical point that one never gets to because it is preempted by the weak first-order transition, and de Gennes, noting this,25 suggested that there ought to be the pre-transitional behavior one would expect from the approach to such a critical point—rapidly increasing correlation lengths and critical slowing down. Indeed, the existence of some sort of meso-scale “pseudo-nematic” ordering is commonly thought of in these terms.16,17

The time-honored way of characterizing the extent of nematic order is to look at the lowest non-vanishing level of alignment with respect to the preferred (director) axis n^. For N molecules whose individual molecular orientation unit vectors are Ω^j (j = 1, …, N), but which have no net head-or-tail preference, the average value of the second Legendre polynomial P2x=123x21 of the components along the director

s0=1Nj=1NP2n^Ω^j
(1.1)

is the standard order parameter.2 Here the brackets refer to an average over all the molecular orientation vectors. This s0 parameter serves as a natural order parameter because it equals 1 when all molecules are perfectly aligned with the director, decreases as the system becomes less ordered, and drops to 0 when the distribution of orientations becomes isotropic.

However, this parameter is obviously of no use in measuring the extent of the order in the isotropic phase: since there is no preferred axis in that phase, there is no director to compare with. In fact, it is not of much computational use in simulations of nematic phases either. Using it as written would require an a priori knowledge of the director, itself an unknown—and a fluctuating one, at that.

The standard alternative employed in computer simulations is to construct an order parameter tensor,2,26

Q=1Nj=1N123Ω^jΩ^j1,
(1.2)

an object with no dependence on any particular axis. When nematic order is present, it is possible to show that the director is an eigenvector of the average 3 × 3 Q matrix, namely, the one with the largest eigenvalue, and that that eigenvalue is, in fact, s0. But simulations, of course, have finite-size-induced fluctuations about that average.27,28 What we want to discuss in this paper is what one can learn from these fluctuations. That is, we want to examine the probability distribution of this Q tensor, p(Q), when the system is not nematically ordered—when its macroscopic limit is to be isotropic. In particular, we would like to discover what finite-size calculations of p(Q) tell us about the nature of the pre-transitional ordering.

There are substantial precedents in the literature for studying finite-size scaling near the isotropic-nematic transition29–31 (and near weakly first order transitions in other systems).32,33 However, unlike most of those studies, we are not interested so much in the macroscopic limit of the transition as we are in features that are only visible with finite numbers of molecules. To that end, we will be making use of a somewhat unusual combination of microscopic statistical mechanics with ideas from random-matrix theory34 and expressing our answers in terms of the extent to which standard random-matrix theory fails to describe fluctuations in liquid-crystalline order.

The remainder of this paper will be organized as follows: Sec. II will review a few salient results and concepts from random-matrix theory and show how the equilibrium p(Q) distribution reverts exactly to the standard random-matrix form in the macroscopic limit of an isotropic phase, regardless of the details of the intermolecular interactions. Section III will then show how similar calculations enable us to derive similarly exact results for equilibrium behavior in nematic phases and for the time-dependence of this order measure in the isotropic case. Section IV will present numerical results derived from computer simulation of a standard Gay-Berne model35 for liquid-crystal-forming fluids, illustrating how well this random-matrix theory works on both sides of the nematic/isotropic transition and showing how and why it begins to fail as that transition is approached from the isotropic side. We conclude in Sec. V with a discussion of the implications of these observations for the intriguing parallels that have been drawn between supercooled liquids and pre-transitional isotropic liquid crystals.11,12,17–23

The sheer impracticality of computing all of the matrix elements in some physically interesting matrices, such as those associated with certain quantum mechanical Hamiltonians, prompted interest in what is usually regarded as the standard random matrix problem: given a matrix whose elements are random variables, what can one say about the distribution of eigenvalues?34,36 Despite the manifestly statistical way in which this question is posed, the results from random matrix theory have found numerous applications to such explicitly deterministic issues as few-body energy-level spacings37 and reaction rates,38 as well as in many-body contexts.39 

The way random matrix theory is normally formulated is to assume that we know nothing about the values of the individual matrix elements except insofar as the symmetry of the systems imposes interrelations between the elements.40 For our Q tensor of Eq. (1.2), a real 3 × 3 matrix whose rows and columns are indexed by the coordinate axes, we know that the matrix has to be symmetric and have zero trace,

Qμ,ν=Qν,μ(μ,ν=x,y,z),TrQ=μ=x,y,zQμ,μ=32Nj=1NΩjx2+Ωjy2+Ωjz21=0.
(2.1)

The probability distribution of Q can therefore be written as a function of the 6 matrix elements

Q1=Qxx,Q2=Qyy,Q3=Qzz,Q4=Qxy,Q5=Qxz,Q6=Qyz
(2.2)

normalized so that

a=16dQapQ=1,
(2.3)
pQδTrQ×some function FQ1,,Q6.
(2.4)

Now suppose, for the purposes of this section, that we are only considering phases that are macroscopically isotropic; ones without any preferred axes. Invariance with respect to the choice of coordinate systems requires that in whatever way p(Q) depends on Q, it must be invariant with respect to orthogonal transformations (i.e., with respect to the replacement,

QO1QO,
(2.5)

with O being an orthogonal matrix).41 The natural candidates for such invariants are traces of Q to arbitrary powers, TrQ,TrQ2,TrQ3,TrQ4,, because such traces are invariant with respect to any similarity transformation (change of basis) such as Eq. (2.5).40 But elementary linear algebra also implies that any TrQn for n > 3 can be expressed in terms of the corresponding traces for n = 1, 2, and 3.42 Thus, given Eq. (2.1), whatever function F(Q) appears in Eq. (2.4) can only involve TrQ2 and TrQ3.43 

Once one has the p(Q) distribution, random matrix theory then provides a simple route to the distribution of the eigenvalues, ρλ1,λ2,λ3. These traces themselves are easily expressed in terms of the three eigenvalues λ1,λ2,λ3. For example,

TrQ=λ1+λ2+λ3, TrQ2=λ12+λ22+λ32.

Inasmuch as the Qa values (a = 1, …, 6) completely determine the eigenvalues, Eqs. (2.3) and (2.4) tell us that the probability distribution of eigenvalues is simply p(Q) multiplied by the Jacobian factor that converts λ1,λ2,λ3Q1,,Q6, the so-called Vandermonde determinant,44 

ρλ1,λ2,λ3λ1λ2λ1λ3λ2λ3×δλ1+λ2+λ3fλ12+λ22+λ32,λ13+λ23+λ33,
(2.6)
λ1>λ2>λ3,

for some function f.

The analysis to this point is exact, independent of strength of intermolecular correlations, so these expressions will be our point of departure for the statistical mechanical development in Sec. II B. However, it is worth noting that this juncture is also where traditional random matrix theory inserts an assumption that we do not wish to make: that (except for the symmetry requirements just mentioned) the matrix elements are governed by independent, identical, zero-centered Gaussian distributions. If one does impose this, apparently unphysical, requirement, the probability distributions of the order-parameter matrix become remarkably simple, as does P(s), the probability distribution of the order parameter s (=λ1),28,40

pQ=342απ5/2δTrQexpαTrQ2,
(2.7a)
α1=25TrQ2,
(2.7b)
ρλ1,λ2,λ3λ1λ2λ1λ3λ2λ3δλ1+λ2+λ3×expαλ12+λ22+λ32,  λ1>λ2>λ3,
Ps  =  sdλ2λ2dλ3ρs,λ2,λ3  =  6απe6αs2192αs2e32αs2,s>0  =  0,               s<0.
(2.7c)

As we now show, even in a fully interacting liquid, Eqs. (2.7a)–(2.7c) become exact in the limit that the number of molecules N → ∞. However, it is the precise dependence of an interacting system’s convergence with increasing system size N that we want to pursue.

It is revealing to look at the microscopic significance of the two symmetry invariants TrQ2 and TrQ3. Using Eq. (1.2) and remembering the identities

Tr1=3, TrΩ^jΩ^j=1, TrΩ^jΩ^k=Ω^jΩ^k

tells us that

TrQ2=32N2j,k=1NP2Ω^jΩ^k,
(2.8)
TrQ3=18N3j,k,l=1N27Ω^jΩ^kΩ^kΩ^lΩ^jΩ^l9Ω^jΩ^k2+Ω^kΩ^l2+Ω^jΩ^l2+6.
(2.9)

In other words, TrQ2 and TrQ3 are two-body and three-body orientational susceptibilities, respectively. It is easy to verify that, for a completely isotropic system, the only two-particle terms [(j, k) with jk] that contribute to TrQ2 are those in which the particle orientations are correlated. Similarly, with TrQ3, one can verify that in the isotropic case, neither terms with three distinct particles nor terms with two distinct particles give nonzero contributions unless those particles are orientationally correlated.

The implication that immediately follows is that these averages ought to scale in well-defined ways with N, the numbers of molecules. If we call Ncorr the typical number of molecules orientationally correlated with a given molecule, then Eqs. (2.8) and (2.9) predict

TrQ21N2N×Ncorr=NcorrN,
(2.10)
TrQ31N3N×Ncorr2=NcorrN2.
(2.11)

In typical molecular liquids, we find Ncorr to be on the order of the size of a solvation shell or two45 so that Ncorr/N1. Hence our expectation going into the calculations is that TrQ2 should be by far the dominant contribution and that its value will ordinarily control the macroscopic limit of p(Q).

The simplest way to verify this supposition is to recognize that the traditional random matrix result (which does depend on just TrQ2)40 is an essentially Gaussian distribution resulting from the operation of the central limit theorem. All of the elements of the Q matrix, Eq. (1.2), are sums over N quantities that fluctuate about 0. So, while it may not be obvious what the ramifications are of the Tr(Q) = 0 requirement, it would not be completely surprising if the elements had Gaussian distributions in the large N limit. It seems eminently reasonable, then, for us to try to make use of the standard methods employed in deriving the central limit theorem.46 

The probability distribution we want can be written as an average over the molecular orientation unit vectors Ω^j (j = 1, …, N),

pQ=pQ1,,Q6=δQ1+Q2+Q3δQ1QxxδQ2Qyy×δQ4QxyδQ5QxzδQ6Qyz,
(2.12)
Qμν=32Nj=1Nqμνj,qμνj=ΩjμΩjν13δμν(μ,ν=x,y,z),
(2.13)

or, using Fourier representations of the delta functions

pQ=δQ1+Q2+Q3a=1a3612πdkaeikaQaGk,
(2.14)
Gk=exp3i2Nj=1Na=1a36kaqaj,
(2.15)

with the conventions that the qaj with a = (1, …, 6) correspond to (xx, yy, zz, xy, xz, yz), respectively, and that k denotes k1,k2,k4,k5,k6.

The generating function G, though, can be expanded in a standard cumulant (fluctuation) expansion, which automatically orders the terms in powers of (1/N),

Gk=exp3i2Njakaqaj98N2ja,kbkakbδqajδqbk+27i48N3ja,kb,lckakbkcδqajδqbkδqcl+,
(2.16)

where the a, b, c indices can take on the values 1, 2, 4, 5, or 6, and expressions such as

δqajqajqaj

denote orientational fluctuations of the jth molecule.

As we show in  Appendix A, in a macroscopically isotropic system, symmetry dictates that all averages over single q’s vanish

qaj=0

and that the only surviving averages over pairs of q’s are simply related to one another. In particular, we prove that the only non-vanishing sums over particles are all proportional to 2-body orientational susceptibility,

1N2j,k=1Nqajqbk=rabTrQ2,
(2.17)
rab(a,b=1,2,4,5,6)=8/135,a=b=1 or 22/45,a=b=4, 5, or 64/135,(a,b)=(1,2) or (2,1)0,otherwise.
(2.18)

Hence, to leading order in (1/N), and in terms of the α constant we introduced in Eq. (2.7b), the generating function is just

Gk=exp124α4k12+4k224k1k2+3k42+3k52+3k62,
(2.19)

which can be substituted back into Eq. (2.14) and integrated,

pQ=342απ5/2δQ1+Q2+Q3FQ,FQ=exp2αQ12+Q22+Q1Q2+Q42+Q52+Q62.
(2.20)

But since Q32=Q1+Q22, we can also write

FQ=expαQ12+Q22+Q32+2Q42+Q52+Q62,
(2.21)

meaning that in the limit of large N, as we expected, the probability distribution of each of the order-tensor’s matrix elements is indeed Gaussian. Moreover, since

TrQ2=μ,ν=x,y,zQμνQνμ=Q12+Q22+Q32+2Q42+Q52+Q62,

we find that the probability distribution of Q in this limit is precisely the traditional random matrix answer we wrote in Eq. (2.7a), with the constant α fixed by Eqs. (2.7b) and (2.8) to be

α=5N31γ, γ=1Nj,k=1NP2Ω^jΩ^k.
(2.22)

The principal significance of Eqs. (2.20)–(2.22) for us is that they tell what to expect for the extent of order in a finite-sized computer simulation of a macroscopically isotropic liquid. The key point is that, based on the same reasoning we used to arrive at Eq. (2.10), we know to look for the constant α in Eq. (2.22) to be of order N. This expectation, in turn, stems from the fact that the γ in Eq. (2.22) is of the order of Ncorr, which we would normally expect to be an N-independent constant. With an ideal (completely orientationally uncorrelated) system, for example, only j = k terms would contribute, so γ = 1 and α = (5N/3).28 

Regardless of the level of intermolecular correlation, if Eqs. (2.20) and (2.21) hold, P(s), the probability distribution of the order parameter s, will be given by Eq. (2.7c) for large N. Since s in this equation only appears in the combination αs2, that means that even though the system has no order in the macroscopic limit, P(s) will not become a delta function at s = 0 until N → ∞. Indeed, the presence of this particular combination implies that the average value of s should scale with system size as

s1/N.
(2.23)

We illustrate the basic scaling behavior in Fig. 1 by considering the aforementioned ideal liquid crystal and comparing the order parameter distribution predicted by Eq. (2.7c) with that obtained in a numerically exact simulation. Consistent with our discussion, even this completely uncorrelated example exhibits a finite amount of order at any finite N because of fluctuations, and the dependence of that order on N is given quite accurately by the traditional random matrix formula.28 

FIG. 1.

Probability distribution of the orientational order parameter s for completely uncorrelated molecules, illustrating the dependence on N, the number of molecules. The plotted points come from molecular dynamics simulations; the solid lines are the predictions of ideal-gas random matrix theory [Eq. (2.7c) with α = 5N/3].

FIG. 1.

Probability distribution of the orientational order parameter s for completely uncorrelated molecules, illustrating the dependence on N, the number of molecules. The plotted points come from molecular dynamics simulations; the solid lines are the predictions of ideal-gas random matrix theory [Eq. (2.7c) with α = 5N/3].

Close modal

Equally importantly, from our derivation, we see that reasonably strong intermolecular correlations are in no way antithetical to traditional random matrix theory. In fact, the only way that these predictions—and their associated finite-size-scaling requirements—can ever fail to hold is for the number of molecules in correlated clusters to start to be comparable to the number of molecules in the system. We show later in the paper that such breakdowns are the hallmark of pre-transitional behavior in liquid crystals.

The basic random matrix framework we have been discussing actually lends itself to a number of useful generalizations. Suppose we are interested in the joint probability distribution of Q tensors at two different times, 0 and t, pQ(0),Q(t). The list of possible invariants enlarges a bit, but remains short, namely,

TrQ(0)=TrQ(t)=0,TrQ(0)Q(0), TrQ(t)Q(t), and TrQ(0)Q(t),
(3.1)

along with the 4 analogous invariants derived from TrQ3. One therefore expects to obtain a random-matrix prediction for pQ(0),Q(t) involving just the three Q2 invariants listed in Eq. (3.1).

In addition, if we define, by analogy with Eq. (2.22), the intensive time-dependent susceptibility

γt=1Nj,k=1NP2Ω^j0Ω^kt
(3.2)

and recognize that [by analogy with Eq. (2.8)]

TrQ(0)Q(t)=32Nγt,TrQ(0)Q(0)=TrQ(t)Q(t)=32Nγ0=32Nγ,
(3.3)

we are again led to expect that this random matrix theory will be exact in the limit that N → ∞ and that the single function γ(t) will control the magnitude of the fluctuations for any finite N.

Applying much the same statistical mechanical approach that we took in Sec. II lets us confirm these predictions. Assume that we remain interested only in macroscopically isotropic phases and only in looking at large N situations. We can generalize Eq. (2.12) by writing the joint probability distribution at times 0 and t as

pQ(0),Q(t)=δQ1(0)+Q2(0)+Q3(0)δQ1(t)+Q2(t)+Q3(t)×FQ(0),Q(t),FQ(0),Q(t)=a=1a36δQa(0)Qa0b=1b36δQb(t)Qbt,
(3.4)

where, for any t,

Qμνt=32Nj=1Nqμνjt,qμνjt=ΩjμtΩjνt13δμν,a,b=1,,6μν=xx,yy,zz,xy,xz,yz.

Using the same Fourier representations of the delta functions,

FQ(0),Q(t)=a=1a3612πdkaeikaQa(0)×b=1b3612πdKbeiKbQb(t)Gk,K,Gk,K=exp3i2Nj=1Na=1a36kaqaj0+b=1b36Kbqbjt,
(3.5)

and performing the same kind of cumulant expansion yields an expression for the generating function. To leading non-vanishing order in (1/N),

Gk,K=exp98N2ja,kbkakbqaj0qbk0+KaKbqajtqbkt+2kaKbqaj0qbkt+.
(3.6)

Once again, symmetry reveals that, in an isotropic system, these averages are all proportional to a 2-body orientational susceptibility [with the same coefficients given in Eq. (2.18)],47 

1N2j,k=1Nqajt1qbkt2=rabTrQt1Qt2,
(3.7)

but where, now, t1 and t2 are two, possibly different, times. The generating function therefore continues to take on a simple form

Gk,K=expγt20N4k1K1+4k2K22k1K22K1k2+3k4K4+3k5K5+3k6K6GkGK

[with Gk and GK referring to the same generating function defined by Eqs. (2.19) and (2.22)]. This expression can be substituted into Eq. (3.5) and integrated to yield a probability distribution depending on just the invariants anticipated in Eq. (3.1),

FQ(0),Q(t)=342απ5Γ5/2texpαΓtSQ(0),Q(t),SQ(0),Q(t)=TrQ(0)Q(0)+TrQ(t)Q(t)2CtTrQ(0)Q(t),Γt=1C2t1
(3.8)

with α given by Eq. (2.22), and C(t) being the normalized collective orientational correlation function,

Ct=γtγ=j,k=1NP2Ω^j0Ω^ktj,k=1NP2Ω^jΩ^k.
(3.9)

When times are long compared to the relaxation time for this function, C(t) → 0 and Γ(t) → 1, so the two-time probability distribution, Eqs. (3.4) and (3.8), reduces, as it must, to the product of the single-time distributions derived in Sec. II,

pQ(0),Q(t)pQ(0)pQ(t),

but, as with single-time distribution functions, the two-time random-matrix result is indeed exact for all times in the large-N limit.

When a liquid loses its macroscopic isotropy, one no longer requires the pQ distribution to be invariant to all possible rotations and exchanges of coordinate axes. In particular, in any system that has macroscopic nematic order, there is a unique axis, the director n^, and a legitimate probability distribution will only be invariant with respect to rotations about that axis.

To put it another way, since a nematic pQ has to distinguish between projections of Q parallel to and perpendicular to the director, n^n^Q and 1n^n^Q, respectively, it cannot be a function of just the TrQn quantities we had been using. Moreover, since the average Q is no longer zero,

Q=s0n^n^s021n^n^,
(3.10)

where s0 is the usual nematic order parameter defined in Eq. (1.1), the invariant ingredients in pQ are most easily constructed based on the fluctuations δQ=QQ instead of the Q tensors themselves. Thus, the tensor that ends up being involved is the most general tensor that has these properties,

P=a1n^n^δQ1n^n^+bn^n^δQ1n^n^+1n^n^δQn^n^+cn^n^δQn^n^,
(3.11)

for some constants a, b, and c. When a = b = c, P is simply proportional to δQ, but in the nematic phase, these three constants will all be different.

Consistent with this observation (as we show in  Appendix B), there are three fundamentally distinct δQ2 invariants in the nematic case,

I1=Trn^n^δQ2,I2=Tr1n^n^δQ2,I3=Trn^n^δQ1n^n^δQ.
(3.12)

Since any positive-definite linear combination of these invariants can be written as TrP2 for some choice of a, b, and c,

TrP2=c2I1+a2I2+2b2I3,
(3.13)

the random-matrix prediction for nematics is of the remarkably simple form

pQδTrQexpα0TrP2,
(3.14a)
α01=25TrP2.
(3.14b)

We can now carry out our statistical mechanics development, as we did before, to show that this form really is the exact large-N result and to arrive at microscopic expressions for the unknown constants.

The setup of the calculation of pQ is identical to what was presented in Sec. II, Eqs. (2.12)–(2.16), but differences appear once we begin to compute the values of the cumulants that enter the generating function, Eq. (2.16). For notational simplicity, suppose we adopt the convention that the director points along the z axis. We find that some of the single-q averages are now finite,

q1j=q2j=13s0, q4j=q5j=q6j=0,

and that some of the degeneracy of the double-q averages is lifted,

1Nj,k=1Nδq1jδq1k=1Nj,k=1Nδq2jδq2k=m1,1Nj,k=1Nδq1jδq2k=m2,1Nj,k=1Nδq4jδq4k=m3,1Nj,k=1Nδq5jδq5k=1Nj,k=1Nδq6jδq6k=m4.
(3.15)

Three of these averages, however, are related by symmetry considerations,

m1m2=2m3,
(3.16)

making it possible to express all four of them in terms of averages of the three invariants given in Eq. (3.12),

m1=N912I1+I2,m2=N932I1I2,m3=N912I1+I2,m4=2N9I3.
(3.17)

Details are given in  Appendix B.

Substituting these results into Eq. (2.16) gives us the nematic version of the generating function to leading order in 1/N,

Gk=exp98Nm1k12+m1k22+2m2k1k2+m3k42+m4k52+m4k62expi2s0k1+k2,

and integrating via Eq. (2.14) yields the desired nematic generalization of Eqs. (2.20) and (2.21),

pQ=2N9π5/2m12m22m3m421/2×δQ1+Q2+Q3FQ,
FQ=exp2Na2Q1+s022+a2Q2+s022+c2Q3s02×  exp4Na2Q42+b2Q52+b2Q62=exp2NTrP2,
(3.18)

in perfect accord with Eqs. (3.11)–(3.14), with the constants a, b, and c now identified as

a2=1912m3=191m1m2,b2=1912m4,c2=19m2m12m22.
(3.19)

In particular, the reader can confirm that Eqs. (3.12)–(3.14), (3.17), and (3.19) imply that

α0=52TrP21=525/4N1=2N,

as required.

Once we have these results, it is also straightforward to arrive at the leading-order-in-N nematic generalization of the eigenvalue and order parameter distributions, Eq. (2.7),

ρλ1,λ2,λ3λ1λ2λ1λ3λ2λ3×δλ1+λ2+λ3exp2Nc2λ1s02+a2λ2+s022+a2λ3+s022
(3.20)
λ1>λ2>λ3,
Ps=sdλ2λ2dλ3ρs,λ2,λ3.

In deriving Eq. (3.20), we simply rewrote Eqs. (3.12)–(3.14) with δQ in its eigenvector basis, making the assumption that the director is an eigenvector for every liquid configuration (which makes the I3 invariant identically zero). This assumption is not, in fact exact, but as we shall see in Sec. IV, the effects of fluctuations in the director direction are numerically infinitesimal in the nematic phase.

Finally, we note that in the isotropic limit, this same equation reverts back precisely to the isotropic answer for the order tensor probability distribution we derived in Sec. II. In that limit, all three invariants contribute, the average order parameter s0 = 0, and the a, b, c coefficients become identical functions of the orientational susceptibility γ defined by Eq. (2.22). The specifics of this limit are also presented in  Appendix B.

The most straightforward way to test these ideas is by simulation of a well-studied model liquid crystal system. Given such a simulation, we can easily calculate the exact distributions of order parameters P(s) for a variety of different system sizes, in both isotropic and nematic phases. We can then make comparisons with the corresponding random matrix predictions by using the same simulations to compute the only parameters that enter those predictions—the average values of the invariants.

All of the simulations presented here are for systems of N molecules interacting via the Gay-Berne potential (N = 256, 500, 864, 1372, and 2048). Gay-Berne potentials are essentially Lennard-Jones representations of the interactions of symmetric-top ellipsoidal molecules, and, as such, are parameterized by the values of the Lennard-Jones diameters σ and well-depths ε for pairs of molecules arranged side-by-side and end-to-end.35 We use one of the standard parameter choices, one whose phase diagram has been extensively mapped,48 

σ||=3σ,σ=σ,εs=5εe,εe=13ε,μ=2,ν=1
(4.1)

(in the usual notation), and limit ourselves to studying a range of dimensionless densities ρσ3 at a fixed dimensionless temperature,

kBT/ε=1.00.

Simulations are carried out via a microcanonical, leapfrog, molecular dynamics algorithm using cubic periodic boundary conditions and largely the same initializing, relaxing, and averaging protocols described in our earlier work.49,50 In particular, our time step was δt=0.001τ0 with τ0=mσ2/ε1/2 (with m being the molecular mass) and averages were computed from data recorded every 10 δt. For the most part, all of our results involving densities ρσ3=0.304 were averaged over 1.5×107 liquid configurations and all of results pertaining to other densities or to uncorrelated systems were averaged over 3×106 liquid configurations. Those instances in which different levels of averaging were employed are explicitly noted in the relevant figure captions.

Orientational ordering information was derived by computing the Q(R) tensor [Eq. (1.2)] at each liquid configuration R we sampled. The matrix was then diagonalized by standard methods and the order parameter for that configuration s(R), identified as the largest eigenvalue. When needed for the nematic case, we also identified the director for that configuration n^R as the corresponding eigenvector.51 

One further practical point is worth mentioning. It is of considerable importance for our subsequent discussion that we can make use of our finite-size scaling ideas to ensure that we are where we think we are on the phase diagram for this model. As we note in Sec. II, for a finite sample, the average order parameter s takes on a non-zero value even in the isotropic phase; the only robust distinction between isotropic and nematic phases is the contrast between the zero and finite values (respectively) that this quantity has in the macroscopic limit. But if finite-size scaling holds, the specific size dependence for isotropic phases, Eq. (2.23) and the comparable equation for nematics,52 allows us to unambiguously characterize the thermodynamic phase (or phases) for any choice of density by extrapolating the average order parameter to infinite system size. Figure 2 shows that this level of scaling does indeed hold quantitatively for all the examples we consider, and that we can safely say that our simulations with reduced densities ρσ30.304 correspond to single-phase macroscopically isotropic liquids, while those with reduced densities ρσ30.33 correspond to single-phase macroscopically nematic liquids.

FIG. 2.

Finite-size scaling behavior of s, the average value of the orientational order parameter s, for a Gay-Berne liquid simulated at three different values of the reduced density ρσ3. N is the number of molecules in each simulation and the points shown all represent averages over 3×106 liquid configurations. The dashes are straight-line fits to the simulation-derived points.

FIG. 2.

Finite-size scaling behavior of s, the average value of the orientational order parameter s, for a Gay-Berne liquid simulated at three different values of the reduced density ρσ3. N is the number of molecules in each simulation and the points shown all represent averages over 3×106 liquid configurations. The dashes are straight-line fits to the simulation-derived points.

Close modal

The N-dependence of the probability distribution of order parameters P(s) for Gay-Berne liquids provides a reasonably robust test of our analysis. Both for a density that corresponds to a macroscopically isotropic state (Fig. 3) and for a density that corresponds to a macroscopically nematic state (Fig. 4), we see that the order-parameter fluctuations are well predicted by random-matrix theory—with an accuracy that clearly increases as N increases. Perhaps the success of a central-limit theorem result is not all that surprising, but our statistical mechanical derivations tell us that to achieve this success, correlations of order-parameter fluctuations about their global averages must be decaying over distances noticeably smaller than the system’s length, even for modest simulation sizes.

FIG. 3.

Probability distribution of the orientational order parameter s for a Gay-Berne liquid simulated at a relatively low density ρ, one that corresponds to a macroscopically isotropic phase. N is the number of molecules in each simulation. The plotted points come from molecular dynamics simulations; the solid lines are the predictions of the random matrix theory presented in Sec. II.

FIG. 3.

Probability distribution of the orientational order parameter s for a Gay-Berne liquid simulated at a relatively low density ρ, one that corresponds to a macroscopically isotropic phase. N is the number of molecules in each simulation. The plotted points come from molecular dynamics simulations; the solid lines are the predictions of the random matrix theory presented in Sec. II.

Close modal
FIG. 4.

Probability distribution of the orientational order parameter s for a Gay-Berne liquid simulated at a relatively high density ρ, one that corresponds to a macroscopically nematic phase. N is the number of molecules in each simulation. The plotted points come from molecular dynamics simulations; the solid lines are the predictions of the random matrix theory presented in Sec. III.

FIG. 4.

Probability distribution of the orientational order parameter s for a Gay-Berne liquid simulated at a relatively high density ρ, one that corresponds to a macroscopically nematic phase. N is the number of molecules in each simulation. The plotted points come from molecular dynamics simulations; the solid lines are the predictions of the random matrix theory presented in Sec. III.

Close modal

In those terms, what is interesting is that as one approaches the isotropic/nematic phase transition from the isotropic side, the story changes. There is a range of intermediate densities within the isotropic regime where random matrix theory seems to fare less well for the very same N values (Fig. 5). The theory still converges to the exact answer as the system size increases, but one evidently has to go to significantly larger sizes to achieve the same level of accuracy (Fig. 6).

FIG. 5.

Probability distribution of the orientational order parameter s for a Gay-Berne liquid simulated over a range of densities ρ, each with N = 1372 molecules. The highest density corresponds to a macroscopically nematic phase; the 3 lower densities correspond to macroscopically isotropic phases. As in Figs. 3 and 4, plotted points are molecular dynamics results, solid lines are random-matrix predictions.

FIG. 5.

Probability distribution of the orientational order parameter s for a Gay-Berne liquid simulated over a range of densities ρ, each with N = 1372 molecules. The highest density corresponds to a macroscopically nematic phase; the 3 lower densities correspond to macroscopically isotropic phases. As in Figs. 3 and 4, plotted points are molecular dynamics results, solid lines are random-matrix predictions.

Close modal
FIG. 6.

Probability distribution of the orientational order parameter s for a Gay-Berne liquid simulated at an intermediate density ρ that still corresponds to a macroscopically isotropic phase. N is the number of molecules in each simulation. Plotted points are molecular dynamics results; solid lines are random-matrix predictions.

FIG. 6.

Probability distribution of the orientational order parameter s for a Gay-Berne liquid simulated at an intermediate density ρ that still corresponds to a macroscopically isotropic phase. N is the number of molecules in each simulation. Plotted points are molecular dynamics results; solid lines are random-matrix predictions.

Close modal

The origin of this phenomenon is nicely captured by a random matrix perspective. According to Eqs. (2.10) and (2.11), random matrix theory’s success in the isotropic phase arises because TrQ2 is the only matrix invariant that contributes significantly for large systems; it scales as N−1, whereas TrQ3 scales as N−2. But those same equations suggest that if the correlation volume ever did become of the order of the system’s volume (if Ncorr ever started to approach N), the scaling would start to fail and the dominance of TrQ2 would no longer be guaranteed. That is exactly what we see in Fig. 7. The N−2 scaling of TrQ3 is maintained only as long as ρσ30.30, and densities ρσ3>0.30 are precisely where we see the random matrix formulas begin to converge slowly.

FIG. 7.

The origin of the failure of finite-size scaling for macroscopically isotropic Gay-Berne liquids at intermediate densities. The two-particle and three-particle isotropic invariants, TrQ2andTrQ3, respectively, are fitted to the indicated power law dependence on N, the number of molecules. The resulting power p is plotted versus reduced density ρσ3 and compared with the 1/N and 1/N2 dependences expected for these invariants (shown as horizontal black lines). The growing loss of dominance of the two-particle invariant at densities ρσ30.30 leads to increasingly prominent failures of finite-size scaling. (The open and filled squares denote the results of alternative fits derived from the different N ranges shown. Each square represents an average over 1.5×107 liquid configurations.)

FIG. 7.

The origin of the failure of finite-size scaling for macroscopically isotropic Gay-Berne liquids at intermediate densities. The two-particle and three-particle isotropic invariants, TrQ2andTrQ3, respectively, are fitted to the indicated power law dependence on N, the number of molecules. The resulting power p is plotted versus reduced density ρσ3 and compared with the 1/N and 1/N2 dependences expected for these invariants (shown as horizontal black lines). The growing loss of dominance of the two-particle invariant at densities ρσ30.30 leads to increasingly prominent failures of finite-size scaling. (The open and filled squares denote the results of alternative fits derived from the different N ranges shown. Each square represents an average over 1.5×107 liquid configurations.)

Close modal

We can look a little more deeply into the intermediate-density region by examining the interplay between the growth of the correlation length and the failure of random-matrix finite-size scaling. Spatially resolved orientational correlations in systems of nonpolar linear molecules can be represented via a pair correlation function such as

Cr=1Nρj,k=1jkNδrrjkP2Ω^jΩ^k,
(4.2)

a function which vanishes at large r for any macroscopic isotropic liquid.20,48 We expect, in particular, that in three dimensions, such spatially resolved correlation functions should asymptotically take on the classical Ornstein-Zernike form

Crer/ξ/r,
(4.3)

with ξ being the correlation length.53 Since a glimpse at the asymptotic behavior of Eq. (4.3) for our liquid crystal system confirms that this equation does represent the decay of correlations in the macroscopically isotropic regime (Fig. 8), we can determine values of ξ from linear fits of lnrCr versus r (Table I).54 

FIG. 8.

Distance dependence of macroscopically isotropic-liquid orientational pair correlation functions and their asymptotic behaviors. The insert shows, on a linear scale, the Cr correlation function versus reduced distance r/σ at an intermediate density, ρσ3=0.304. The main panel shows the same kind of correlation function plotted as lnrCr over a range of densities. For both, we compare the simulated correlation functions to the results of a fit to Eq. (4.3) (dashed lines). All of these correlation functions are derived from simulations with N = 2048 molecules and averaged over 5×104 liquid configurations.

FIG. 8.

Distance dependence of macroscopically isotropic-liquid orientational pair correlation functions and their asymptotic behaviors. The insert shows, on a linear scale, the Cr correlation function versus reduced distance r/σ at an intermediate density, ρσ3=0.304. The main panel shows the same kind of correlation function plotted as lnrCr over a range of densities. For both, we compare the simulated correlation functions to the results of a fit to Eq. (4.3) (dashed lines). All of these correlation functions are derived from simulations with N = 2048 molecules and averaged over 5×104 liquid configurations.

Close modal
TABLE I.

Single-molecule and collective orientational properties of isotropic liquid crystal states.a

Single-molecule propertiesbCollective propertiesc
ρσ3γ1st-shellγτ/τ0ξ/σNcorr
0.29 2.48 7.19 19.32 1.98 2.24 
0.30 2.70 14.99 45.63 3.15 9.38 
0.304 2.90 28.47 104.59 12.88 650 
Single-molecule propertiesbCollective propertiesc
ρσ3γ1st-shellγτ/τ0ξ/σNcorr
0.29 2.48 7.19 19.32 1.98 2.24 
0.30 2.70 14.99 45.63 3.15 9.38 
0.304 2.90 28.47 104.59 12.88 650 
a

Orientational susceptibilities (γ), correlation lengths (ξ), and relaxation times (τ) derived from simulating a single temperature (kBT/ε = 1.00) Gay-Berne liquid at different densities ρ, all of which would lead to an isotropic phase in the macroscopic limit. Susceptibilities and correlation lengths come from simulations with N = 2048 molecules; relaxation times are from N = 1372 simulations. Averages are over 1.5×107 configurations for susceptibilities, and over 1.2×107 and 5×104 configurations for the correlation functions used to define relaxation times and correlation lengths, respectively. The basic Gay-Berne units of energy, distance, and time (ε, σ, and τ0) are described in the text.

b

The first-shell susceptibility is given by Eq. (2.22) with molecular pairs (j, k) restricted to those within the first-shell distance determined by the first peak of the center-of-mass radial distribution function.

c

The collective orientational susceptibility γ is defined by Eq. (2.22). The relaxation times τ and correlation lengths ξ are determined from log-linear fits to Eq. (3.9) and to Eqs. (4.2)–(4.3), respectively. The effective number of correlated molecules Ncorr=ρξ3.

With these numbers in hand, it is then interesting to look quantitatively at the deviations from random-matrix behavior. What random matrix theory is essentially predicting is that macroscopic orientational order ought to be governed by the central limit theorem and therefore must have a Gaussian distribution in the large N limit. We can test this basic premise by measuring a non-Gaussian parameter, just as one does when investigating non-Gaussian features of supercooled liquids.55,56 If the basic random matrix equation, Eq. (2.7a), holds, then we can define a generating function,

Zα=dQδTrQexpαTrQ2=43π2α5/2,

and use it to prove that in the large-N limit

TrQ22=1Z2Zα2=75lnZα2=75TrQ22.
(4.4a)

Hence, the dimensionless non-Gaussian parameter

α2=TrQ2275TrQ22TrQ22
(4.4b)

should vanish as the system size increases.

As is evident from Fig. 9, before the intermediate density regime is reached, this α2 parameter does indeed decrease monotonically with increasing N. That trend is precisely what we would have expected from the increasing accuracy of random matrix theory we saw in Fig. 3. However, for ρσ3=0.304, a density in the intermediate regime, α2 has a pronounced maximum around N = 1000. It is apparently only for larger N values that an eventual decay toward 0 begins to occur and random matrix theory begins to become accurate here.

FIG. 9.

The dependence of the value of the static non-Gaussian parameter α2 on N, the number of molecules (256 ≤ N ≤ 2048) in macroscopically isotropic Gay-Berne liquids. The relatively low-density liquid (ρσ3=0.29) monotonically approaches the fully Gaussian, random-matrix, behavior (α2 = 0) as the system size increases. However the intermediate density liquid (ρσ3=0.304) does not begin its decay toward random-matrix behavior until after the system size surpasses a critical size. Shown for comparison, as dashed vertical lines, are N=Ncorr, the number of molecules in a correlation volume ξ3, at this intermediate density (with the necessary correlation length ξ determined from the asymptotic decay of the relevant distance-dependent orientational correlation function shown in Fig. 8) and N=2Ncorr. All non-Gaussian parameters reported represent averages over 1.5×107 liquid configurations.

FIG. 9.

The dependence of the value of the static non-Gaussian parameter α2 on N, the number of molecules (256 ≤ N ≤ 2048) in macroscopically isotropic Gay-Berne liquids. The relatively low-density liquid (ρσ3=0.29) monotonically approaches the fully Gaussian, random-matrix, behavior (α2 = 0) as the system size increases. However the intermediate density liquid (ρσ3=0.304) does not begin its decay toward random-matrix behavior until after the system size surpasses a critical size. Shown for comparison, as dashed vertical lines, are N=Ncorr, the number of molecules in a correlation volume ξ3, at this intermediate density (with the necessary correlation length ξ determined from the asymptotic decay of the relevant distance-dependent orientational correlation function shown in Fig. 8) and N=2Ncorr. All non-Gaussian parameters reported represent averages over 1.5×107 liquid configurations.

Close modal

We suggest that the key to understanding this dichotomy is to realize that this threshold for Gaussian behavior is of the order of Ncorr=ρξ3, the number of molecules in a correlation volume (Table I). Figure 9 shows that this claim is correct for intermediate densities, but this observation is also consistent with what we see for lower densities, where Ncorr is smaller than the smallest simulation attempted, resulting in the absence of a visible threshold. The basic idea is that for N>Ncorr, the system is made up of multiple pseudo-nematic domains with lengths on the order of ξ. If there are enough of these domains, the strong correlations within the domains are irrelevant to the global order because the relative scalings of the two- and three-particle matrix invariants, Eqs. (2.10) and (2.11),

TrQ2TrQ3NNcorr=the number of domains,

will always cause the two-particle invariant to dominate, leading to a Gaussian distribution of orientational order. However, if NNcorr, then the system acts as a single domain and the probability distribution of order parameters for a putatively isotropic system can extend well into the nematic range (Fig. 6).

The existence of pseudo-nematic domains has profound consequences for the dynamics of liquid crystals approaching the isotropic side of the isotropic/nematic transition. While the single-molecule analog of the C(t) collective reorientational time correlation function introduced in Eq. (3.9),

C1t=1Nj=1NP2Ω^j0Ω^jt,
(4.5)

shows no significant trends in relaxation rates as the transition is approached, the collective function exhibits a dramatic slowing down with increasing density (Fig. 10 and Table I).19,57

FIG. 10.

Orientational relaxation of macroscopically isotropic Gay-Berne liquids displayed as a function of reduced time t/τ0. Plotted in the main panel are the normalized collective orientational time correlation functions, Eq. (3.9), which show a pronounced growth in relaxation times as the density ρ is increased. By way of comparison, the inset shows that the corresponding single-molecule orientational time correlation functions, Eq. (4.5), relax much more quickly, with far less dependence on thermodynamic conditions. All of these curves are derived from simulations with N = 2048 molecules. The collective and single-molecule correlation functions are averaged over 1.2×107 and 7000 liquid configurations, respectively.

FIG. 10.

Orientational relaxation of macroscopically isotropic Gay-Berne liquids displayed as a function of reduced time t/τ0. Plotted in the main panel are the normalized collective orientational time correlation functions, Eq. (3.9), which show a pronounced growth in relaxation times as the density ρ is increased. By way of comparison, the inset shows that the corresponding single-molecule orientational time correlation functions, Eq. (4.5), relax much more quickly, with far less dependence on thermodynamic conditions. All of these curves are derived from simulations with N = 2048 molecules. The collective and single-molecule correlation functions are averaged over 1.2×107 and 7000 liquid configurations, respectively.

Close modal

That this effect really arises from the presence of growing pseudo-nematic domains is, again, something that random matrix theory helps us verify. We have been suggesting that a failure of random-matrix finite-size scaling in the equilibrium order parameter distribution is indicative of growing pseudo-nematic domains. But, a scaling failure also appears in the dynamics. The single-molecule correlation function Csinglet certainly looks as if it should not depend on N and, making use of the same kind of analysis we did in Sec. II, leads us to expect that both γ(t) and γ entering into the collective version, Ct=γt/γ, Eq. (3.9), should be N-independent quantities proportional to Ncorr. Yet, a closer look at the collective correlation function, Fig. 11, shows that, near the transition, it does have an N dependence, with the multiple domain N = 2048 case relaxing noticeably faster than the single-domain N = 256, 500, and 864 cases for ρσ3=0.304.58 

FIG. 11.

Tests of finite-size scaling for the collective orientational time correlation functions of macroscopically isotropic Gay-Berne liquids. The top and bottom panels (respectively) examine the dependence of these correlation functions on the number of molecules, N, for a relatively low density, ρσ3=0.29, and an intermediate density, ρσ3=0.304. All of these correlation functions are averaged over 1.2×107 liquid configurations.

FIG. 11.

Tests of finite-size scaling for the collective orientational time correlation functions of macroscopically isotropic Gay-Berne liquids. The top and bottom panels (respectively) examine the dependence of these correlation functions on the number of molecules, N, for a relatively low density, ρσ3=0.29, and an intermediate density, ρσ3=0.304. All of these correlation functions are averaged over 1.2×107 liquid configurations.

Close modal

As with the time-independent distribution functions, the deviations from the expected finite-size scaling seen here are indicative of deviations from the Gaussian behavior predicted by random-matrix theory. It is therefore of some use to create a dynamical non-Gaussian parameter analogous to the static version of Eq. (4.4).22,23,55,56 Inasmuch as the random-matrix predictions for the two-time probability distribution pQ(0),Q(t) have the form shown in Eqs. (3.4) and (3.8), we can use the generating function,

ZA,B=dQ(0)dQ(t)δTrQ(0)δTrQ(t)×expαSA,B,
SA,B=ATrQ(0)2+TrQ(t)22BTrQ(0)Q(t),

to compute the corresponding random-matrix-theory (RMT) averages at the relevant RMT values: A=Γ(t),B=C(t)Γ(t). Since ZA,Bα5A2B25/2,

TrQ(0)Q(t)2=12αlnZA,BBRMT2=TrQ22C2t,
TrQ(0)Q(t)2=14α2Z2ZA,BB2RMT=TrQ2215+65C2t.

Because the collective time correlation function C(t) approaches 0 as t → ∞, we can express these quantities in terms of 2-point and 4-point time correlation functions, respectively,

C2t=Ct=TrQ(0)Q(t)TrQ2,
C4(t)=TrQ(0)Q(t)2TrQ(0)Q()2TrQ22,

which also vanish at long times. (Here, the ∞ superscript refers to the limit t → ∞.) Hence we can write the time-dependent analog of Eq. (4.4b) by defining a dynamical non-Gaussian parameter α2t,

α2t=C4t65C22t.
(4.6)

Since C(0) = 1, it is easy to confirm that α20=α2, the static non-Gaussian parameter of Eq. (4.4b), and since α2=0, the decay of α2t must be capturing the relaxation dynamics of the pseudo-nematic domains–whose existence on length scales comparable to the size of our system leads (in the view of our random-matrix theory) to deviations of the two-time probability distribution, pQ(0),Q(t), from a Gaussian form. We shall consider the quantitative behavior of this α2t in future work, but we note here the prominent role of “two-molecules-each-considered-at-two-times” (j = k, l = m) kinds of terms in our 4-point correlation function

TrQ(0)Q(t)2=32N22j,k,l,m=1NP2Ω^j0Ω^kt×P2Ω^l0Ω^mt,

precisely the kind of correlations crucial to the χ4 statistic used to understand supercooled liquids.59,60

What is surprising about having pseudo-nematic domains appear in an orientational disordered phase is not so much that such domains exist but how dramatically they grow in as the isotropic/nematic transition is approached. Both correlation lengths and collective relaxation times seem to follow the Landau-de Gennes prediction that they act as if they were approaching a critical point, with a divergence at some effective critical density ρc,

ξ1/ρcρ1/2,τ1/ρcρ,
(4.7)

that would lie within the nematic region if it were ever reached.2,25,61,62

It is actually straightforward to show that the Gay-Berne model we have been using exhibits just this phenomenology. We could extract correlation lengths from a lengthy series of fits of the kind displayed in Fig. 8, but since we only need relative values, it is easier to note that an integral over the orientational pair correlation function in Eq. (4.2) yields the orientational susceptibility γ defined by Eq. (2.22),

γ=ρdrCr,

which itself is simply related to the correlation length. Integrating the asymptotic form, Eq. (4.3), for example, would imply that γ=4πξ2 and hence that we should expect

γ1/ρcρ.
(4.8)

We show in Fig. 12 that our susceptibilities and collective relaxation times follow Eqs. (4.7) and (4.8) extremely well.

FIG. 12.

The divergences of orientational correlation lengths and collective orientational relaxation times with increasing density ρ. The collective orientational relaxation times, τ, are extracted from the asymptotic behavior of collective orientational correlation functions and plotted as relaxation rates τ−1 (rightmost points). The correlation lengths ξ are not shown directly, but since the collective orientational susceptibilities γ ∼ ξ2, the behavior of ξ−2 is captured by the reciprocal susceptibilities, γ−1 (leftmost points). Independent straight lines fitted through both sets of points (dashed lines) demonstrate the existence of Landau-de Gennes power law divergences at very similar critical densities (ρσ3=0.307 and 0.309, respectively). Data for this figure were collected from simulations with N = 1372 molecules. Averaging was carried out using 1.2×107 and 1.5×107 liquid configuration for τ and γ, respectively.

FIG. 12.

The divergences of orientational correlation lengths and collective orientational relaxation times with increasing density ρ. The collective orientational relaxation times, τ, are extracted from the asymptotic behavior of collective orientational correlation functions and plotted as relaxation rates τ−1 (rightmost points). The correlation lengths ξ are not shown directly, but since the collective orientational susceptibilities γ ∼ ξ2, the behavior of ξ−2 is captured by the reciprocal susceptibilities, γ−1 (leftmost points). Independent straight lines fitted through both sets of points (dashed lines) demonstrate the existence of Landau-de Gennes power law divergences at very similar critical densities (ρσ3=0.307 and 0.309, respectively). Data for this figure were collected from simulations with N = 1372 molecules. Averaging was carried out using 1.2×107 and 1.5×107 liquid configuration for τ and γ, respectively.

Close modal

What this paper has been concerned with are the fluctuations that occur in global measures of liquid-crystalline orientational order, and since the fractional contribution of those fluctuations has to vanish in the limit that the number of molecules N becomes infinite, the paper has essentially been about how the study of finite-size samples (and finite-size domains within a sample) reveals those fluctuations.

From that perspective, one of the principal lessons we learned here is that we can add the description of fluctuations in liquid crystalline order to the list of strongly interacting many-body problems that can be thought of in random matrix terms. Although random-matrix representations are conventionally regarded as phenomenological, we have shown how one can use straightforward statistical mechanics to explain why the equilibrium probability distributions for liquid crystalline order should follow random-matrix predictions exactly for large N, and should do so in both the isotropic and nematic phases. In fact, we have shown that even the probability distribution for the time evolution of the order has to obey the dictates of random matrix theory in this limit. The fundamental reason for this agreement, in each case, is that two-body correlations determine the statistics of the fluctuations in the large N limit. The end result is that random-matrix answers (which are little more than the various different Gaussian distributions one would have expected from the operation of the central limit theorem in each situation) nicely capture thermodynamic limit behavior.

The basic validity of this analysis was demonstrated by our observation that random-matrix theory is usually quite accurate in forecasting the N-dependence of simulated orientational-order probability distributions. When the correlation lengths of orientational fluctuations are short, as they are in both the nematic phase and the isotropic phase far from the isotropic/nematic phase boundary, the widths of these distributions display precisely the finite-size scaling expected from phenomenological random matrix theory. However, as the correlation lengths in putatively isotropic liquids start to approach values on the order of the sample size—that is, when there start to be identifiable pseudo-nematic domains in the system—this finite-size scaling begins to fail and the usefulness of understanding the statistical mechanical underpinnings of random-matrix theory becomes evident.

The random matrix framework is still helpful to us in that it tells us that no matter how long-ranged the correlations become, one never needs to study more than the two- and three-body invariants to understand the system’s developing order. But statistical mechanics is what tells us how the growing importance of the three-body invariant affects that order: how the Gaussian statistics is not recovered until the sample size becomes significantly greater than the domain sizes and how the system’s orientational dynamics begin to slow dramatically as the isotropic/nematic (I/N) transition is approached.

The behavior of the dynamics in the pre-transitional region of liquid crystals is especially intriguing in view of the analogies that have been made with supercooled liquid dynamics.11,12,17–23 It has been pointed out that the asymptotic relaxation times seen in optical Kerr measurements on isotropic liquid crystals10–13,15,21 lengthen rapidly as the I/N transition is approached, calling to mind the rapid divergences seen in the viscosities of supercooled liquid en route to their glass transitions. But the observed liquid crystal orientational relaxation times diverge in just the way that Eq. (4.7) (or more precisely, its temperature dependent equivalent) would have predicted, rather than in the exponentially more explosive Vogel-Tammann-Fulcher fashion that characterizes supercooled liquids.59 Moreover, unlike supercooled liquids, which show pronounced slowing down of both their single-molecule, Eq. (4.5), and collective, Eq. (3.9), orientational relaxation times,63,64 with liquid crystals it is only the collective orientational relaxations that slow noticeably. Optical Kerr measurements, in particular, report on collective polarizability time correlation functions. Indeed, for linear molecules examined for times long enough that correlations in the interaction-induced polarizability become irrelevant, such polarizability correlation functions revert precisely to the very collective P2 orientational correlation function we have been using.65 

These observations do not diminish the significance of the similarities between pre-transitional liquid crystals and supercooled liquids, but they do suggest that the essential similarity may lie in the heterogeneities in both systems rather than in any parallel behaviors of the individual degrees of freedom.23,66 We have shown that non-Gaussian fluctuations in liquid crystalline order, for example, are a direct reflection of the existence of a finite number of domains. More intriguingly, any non-Gaussian character to the time dependence of that order reveals the χ4 statistic diagnostic of incipient glassy behavior,59,60 paving the way, perhaps, to future spectroscopic avenues for probing the fundamental origins of the intriguing parallels with glassy systems.

See supplementary material for details concerning the molecular dynamics protocols used to ensure that we had appropriately relaxed liquid configurations.

This work was supported by the National Science Foundation under Grant No. CHE-1565540. The original molecular dynamics code used in this work was created by Layne Frechette. We thank him and Vale Cofer-Shabica for valuable discussions on a variety of computational matters.

Both the developments in Sec. II and in Sec. III A require that we evaluate averages over the single-molecule orientational tensor elements qμνj [introduced in Eq. (2.13)] in isotropic systems (j=1,,N;μ,ν=x,y,z). However, the invariance of isotropic phases with respect to arbitrary rotations and reflections of the laboratory axes means that most of the averages vanish and that the others are degenerate.

In particular, since each orientational vector Ω^j is a unit vector, μΩjμ2=1, and since averages must be invariant with respect to the choice of any single axis μ and with respect to any reflection of the form ΩjμΩjμ, in an isotropic system, we have

ΩjμΩjν=13δμν.
(A1)

Hence, all single-molecule averages must vanish,

qμνj=0.
(A2)

Similarly, averages over the product of two tensor elements qμνjt1qγδkt2 (involving what may be two different molecules j and k at two different times t1andt2) can be evaluated using the isotropic-system identities

3Ωjμ2t1Ωkμ2t2+6Ωjμ2t1Ωkν2t2μν=1,
(A3)
Ωjμ2t1Ωkμ2t2Ωjμ2t1Ωkν2t2μν=2Ωjμt1Ωjνt1Ωkμt2Ωkνt2μν.
(A4)

Equation (A3) follows from the unit-vector property μ,νΩjμ2t1Ωkν2t2=1, and Eq. (A4) can be derived by noting that averages must be invariant to 45° rotations around any third Cartesian axis γμ,ν,

ΩiμtΩiνt12Ωiμt+ΩiνtΩiμtΩiνt,

for any i and t.

Inasmuch as we can take an axis μ to point in any arbitrary direction Ω^,

qμμjt1qμμkt2=49P2Ω^jt1Ω^P2Ω^kt2Ω^=4914πdΩ^P2Ω^jt1Ω^×P2Ω^kt2Ω^=445P2Ω^jt1Ω^kt2.
(A5)

So, we find that, using Eqs. (A1), (A3), and (A4),

qμμjt1qννkt2μν=245P2Ω^jt1Ω^kt2,
(A6)
qμνjt1qμνkt2μν=115P2Ω^jt1Ω^kt2
(A7)

and that all other averages of the form qμνjt1qγδkt2 vanish.

If these results are viewed in conjunction with Eq. (2.8), it is evident that Eqs. (A5)–(A7) are precisely the result expressed in Eqs. (2.17) and (2.18).

The expressions derived in  Appendix A for averages over the single-molecule orientational tensor elements qμνj [Eq. (2.13)] can be generalized to nematic systems. The existence of a unique director axis in uniaxial nematics means that some of the rotational symmetry that characterizes isotropic phases is lost, but expressions governing nematics must still be invariant with respect to rotations about the director axis, as well as to any reflection of the form ΩjμΩjμ. Suppose we consider the consequences of these residual symmetries, calling μ = z the director axis.

As before, the components of each orientational vector Ω^j satisfy the normalization condition, Ωjx2+Ωjy2+Ωjz2=1. Thus, identifying the average order parameter as s0, Eq. (1.1), and using the rotational symmetry in the xy plane tells us that

Ωjz2=132s0+1,Ωjx2=Ωjy2=131s0,

so that

qzzj=23s0, qxxj=qyyj=13s0.
(B1a)

All other single-molecule averages must vanish

qμνjμν=0.
(B1b)

Averages over pairs of tensor fluctuations

δqμνj=qμνjqμνj=ΩjμΩjνδμνΩjμ2

proceed similarly to those in  Appendix A as well, but as with the single-molecule averages, in nematics the (x, y, z) degeneracy is split.

With the isotropic case of random matrix theory, all of these pair averages, Eq. (2.17), and thus the full probability distribution of the Q tensor, depend on a single parameter TrQ2 or, equivalently, on a single susceptibility, γ [Eqs. (2.8) and (2.22)]. The loss of degeneracy in nematics means that there are a larger number of parameters to specify, but that number is still limited by symmetries. The continued presence of reflection symmetries, for example, still means that the only potentially distinct, non-vanishing, averages of the form δqajδqbk (with, what may be, two different molecules j and k and with a, b = xx, yy, xy, xz, yz = 1, 2, 4, 5, 6) will have each x, y, z index appearing an even number of times. That is, (a, b) = (1, 1), (2, 2), (1, 2), (2, 1), (4, 4), (5, 5), or (6, 6). Moreover, as before, the invariance with respect to 90° rotations about the director (z) axis (or equivalently, exchanges xy) makes (1, 1) = (2, 2), (1, 2) = (2, 1), and (5, 5) = (6, 6).

These considerations leave us with four potentially independent parameters, the four susceptibilities m1, m2, m3, m4 defined in Eq. (3.15). In terms of the particular coordinate system we are using,

m1  =  1Nj,k=1NΩjx2Ωkx2Ωjx22,m2  =  1Nj,k=1NΩjx2Ωky2Ωjx22,m3  =  1Nj,k=1NΩjxΩjyΩkxΩky,m4  =  1Nj,k=1NΩjxΩjzΩkxΩkz.
(B2)

However there is one other symmetry we can invoke: the special case of Eq. (A4) requiring that averages must be invariant to 45° rotations around the director (z) axis,

Ωjx2Ωkx2Ωjx2Ωky2=2ΩjxΩjyΩkxΩky.
(B3)

Consequently, we know that the first three of these parameters are related to one another [as noted in Eq. (3.16)],

m1m2=2m3,
(B4)

leaving us with only three independent parameters needed to describe orientational order at the two-body level.

An alternate but equivalent issue is what the orientational invariants are for a nematic system though second order in fluctuations of the order parameter tensor. The requisite quantities now have to be invariant with respect to rotation about the director n^ but (unlike the case with isotropic systems) should not be invariant with respect to rotations about other axes. Thus, the lone orientational invariant for isotropic systems, TrQ2, for example, is no longer appropriate.

A formal way to implement these requirements is to note that on rotating the entire system by an infinitesimal angle about the director n^, the orientation of the jth molecule transforms Ω^jΩ^j+dφn^×Ω^j so that fluctuations in the orientational order parameter tensor,

δQ=32Nj=1Nδqj, δqj=Ω^jΩ^jΩ^jΩ^j,
Ω^jΩ^j=2s0+13n^n^+1s031n^n^,

themselves transform as

δQδQ=32Nj=1Nδqj,
δqj=δqj+dφpj, pj=n^×Ω^jΩ^j+Ω^jn^×Ω^j.

Using this observation, it is straightforward to show that the three quantities introduced in Eq. (3.12) are invariant to infinitesimal rotations about the director n^,

I1=Trn^n^δQ2=Trn^n^δQ2,I2=Tr1n^n^δQ2=Tr1n^n^δQ2,I3=Trn^n^δQ1n^n^δQ=Trn^n^δQ1n^n^δQ.

One can also show, by constructing the analogous transformation formulas for rotations about arbitrary axes, that these expressions are not invariant to infinitesimal rotations about other axes—meaning that they are indeed invariants specific to nematics. Furthermore, since

δQ=n^n^δQ+1n^n^δQ,

any two-body invariant that has components taken from δQ2 can be represented as a linear combination of I1,I2, andI3, Eq. (3.13). Hence these three quantities are the only two-body invariants for nematics. For example,

TrδQ2=I1+I2+2I3.
(B5)

It is easy to see that, within the n^=z^ coordinate system used here,

I1=δQzz2, I2=δQxx2+δQyy2+2δQxy2, I3=δQxz2+δQyz2,
(B6)

telling us, with the aid of Eqs. (2.13) and the Ω^ unit-vector normalization condition, that

I1=94N2j,k=1Nδq3jδq3k=94N2j,k=1Nδq1j+δq2jδq1k+δq2k,I2=94N2j,k=1Nδq1jδq1k+δq2jδq2k+2δq4jδq4k,I3=94N2j,k=1Nδq5jδq5k+δq6jδq6k.
(B7)

The corresponding averages of these invariants can therefore be expressed in terms of the m1,m2,m3,m4 susceptibilities of Eq. (3.15). On invoking the xy invariance and Eq. (B4), we obtain

I1=92Nm1+m2,I2=92Nm1+m3=94N3m1m2,I3=92Nm4,
(B8)

which can be rearranged, if desired, into the results given in Eq. (3.17).

To conclude this appendix, we consider the isotropic limit of all these expressions. In order to have the nematic distribution, Eq. (3.14), revert to the isotropic distribution, Eq. (2.7), we require

α0Tr P2=αTr Q2.

But since Q=δQ in the isotropic case, Eqs. (2.22), (3.13), and (B5) imply that this condition is satisfied if

a2=b2=c2=α/α0=5/6γ.
(B9)

This result, in turn, has implications for the values of the m susceptibilities. If Eq. (3.19) is to hold with the a, b, and c coefficients identical to one another,

m1=2m2, m3=m4,

and thus, in the isotropic limit,

m1=4γ/45, m2=2γ/45, m3=m4=γ/15,NI1=γ/5, NI2=7γ/10, NI3=3γ/10.
(B10)
1.
P. G.
de Gennes
and
J.
Prost
, in
The Physics of Liquid Crystals
, 2nd ed. (
Clarendon
,
Oxford
,
1993
), Chap. 1.
2.
P. G.
de Gennes
and
J.
Prost
, in
The Physics of Liquid Crystals
, 2nd ed. (
Clarendon
,
Oxford
,
1993
), Chap. 2.
3.
J. D.
Litster
and
T. W.
Stinson
 III
,
J. Appl. Phys.
41
,
996
997
(
1970
).
4.
T. W.
Stinson
 III
and
J. D.
Litster
,
Phys. Rev. Lett.
25
,
503
506
(
1970
).
5.
T. W.
Stinson
 III
and
J. D.
Litster
,
Phys. Rev. Lett.
30
,
688
692
(
1973
).
6.
T.
Ueno
,
K.
Sakai
, and
K.
Takagi
,
Phys. Rev. E
54
,
6457
6461
(
1996
).
7.
J. J.
Krich
,
M. B.
Romanowsky
, and
P. J.
Collins
,
Phys. Rev. E
71
,
051712
(
2005
).
8.
Z.
Zhang
,
F.
Wang
, and
Y.
Han
,
Phys. Rev. Lett.
107
,
065702
(
2011
).
9.
C. K.
Mishra
,
K. H.
Nagamanasa
,
R.
Ganapathy
,
A. K.
Sood
, and
S.
Gokhale
,
Proc. Natl. Acad. Sci. U. S. A.
111
,
15362
15367
(
2014
).
10.
H.
Cang
,
J.
Li
, and
M. D.
Fayer
,
Chem. Phys. Lett.
366
,
82
87
(
2002
).
11.
H.
Cang
,
J.
Li
,
V. N.
Novikov
, and
M. D.
Fayer
,
J. Chem. Phys.
118
,
9303
(
2003
).
12.
J.
Li
,
H.
Cang
,
H. C.
Andersen
, and
M. D.
Fayer
,
J. Chem. Phys.
124
,
014902
(
2006
).
13.
J.
Li
,
I.
Wang
, and
M. D.
Fayer
,
J. Chem. Phys.
124
,
044906
(
2006
).
14.
K. P.
Sokolowsky
and
M. D.
Fayer
,
J. Phys. Chem. B
117
,
15060
15071
(
2013
).
15.
K. P.
Sokolowsky
,
H. E.
Bailey
, and
M. D.
Fayer
,
J. Phys. Chem. B
118
,
7856
7868
(
2014
).
16.
K. P.
Sokolowsky
,
H. E.
Bailey
, and
M. D.
Fayer
,
J. Chem. Phys.
141
,
194502
(
2014
).
17.
K. P.
Sokolowsky
,
H. E.
Bailey
,
D. J.
Hoffman
,
H. C.
Andersen
, and
M. D.
Fayer
,
J. Phys. Chem. B
120
,
7003
7015
(
2016
).
18.
D.
Chakrabarti
and
B.
Bagchi
,
Adv. Chem. Phys.
141
,
249
(
2009
).
19.
P. P.
Jose
and
B.
Bagchi
,
J. Chem. Phys.
120
,
11256
(
2004
).
20.
P. P.
Jose
and
B.
Bagchi
,
J. Chem. Phys.
121
,
6978
(
2004
).
21.
D.
Chakrabarti
and
B.
Bagchi
,
J. Chem. Phys.
126
,
204906
(
2007
).
22.
B.
Jana
and
B.
Bagchi
,
J. Chem. Sci.
119
,
343
(
2007
).
23.
B.
Jana
,
D.
Chakrabarti
, and
B.
Bagchi
,
Phys. Rev. E
76
,
011712
(
2007
).
24.
E. F.
Gramsbergen
,
L.
Longa
, and
W. H.
de Jeu
,
Phys. Rep.
135
,
195
(
1986
).
25.
P. G.
de Gennes
,
Phys. Lett. A
30
,
454
(
1969
);
P. G.
de Gennes
,
Mol. Cryst. Liq. Cryst.
12
,
193
(
1971
).
26.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Clarendon
,
Oxford
,
1987
), Sec. 11.5.
27.
R.
Eppenga
and
D.
Frenkel
,
Mol. Phys.
52
,
1303
1334
(
1984
).
28.
T. P.
Doerr
,
D.
Herman
,
H.
Mathur
, and
P. L.
Taylor
,
Europhys. Lett.
59
,
398
402
(
2002
).
29.
M.
Greschek
and
M.
Schoen
,
Phys. Rev. E
83
,
011704
(
2011
).
30.
N. V.
Priezjev
and
R. A.
Pelcovits
,
Phys. Rev. E
63
,
062702
(
2001
).
31.
H.
Weber
,
W.
Paul
, and
K.
Binder
,
Phys. Rev. E
59
,
2168
(
1999
).
32.
P.
Peczak
and
D. P.
Landau
,
Phys. Rev. B
39
,
11932
(
1989
).
33.
K.
Vollmayr
,
J. D.
Reger
,
M.
Scheucher
, and
K.
Binder
,
Z. Phys. B: Condens. Matter
91
,
113
(
1993
).
34.
M. L.
Mehta
,
Random Matrices
, 2nd ed. (
Academic
,
Boston
,
1991
).
35.
J. G.
Gay
and
B. J.
Berne
,
J. Chem. Phys.
74
,
3316
(
1981
).
36.
E. P.
Wigner
,
SIAM Rev.
9
,
1
(
1967
).
37.
E. B.
Stechel
and
E. J.
Heller
,
Annu. Rev. Phys. Chem.
35
,
563
(
1984
).
38.
U.
Peskin
,
W. H.
Miller
, and
H.
Reisler
,
J. Chem. Phys.
102
,
8874
(
1995
).
39.
T.
Halpin-Healy
,
Phys. Rev. Lett.
109
,
170602
(
2012
).
40.
M. L.
Mehta
, in
Random Matrices
, 2nd ed. (
Academic
,
Boston
,
1991
), Chap. 2.
41.

Rotations about arbitrary axes, reflections through arbitrary planes, and permutations of axes can all be represented by orthogonal matrices.

42.
The Cayley-Hamilton theorem guarantees that any p×p matrix satisfies the same characteristic equation as its eigenvalues (a pth order polynomial equation), so any Qn (n > 3) can be expressed as a linear combination of Q,Q2,and Q3. See, for example,
K.
Hoffman
and
R.
Kunze
,
Linear Algebra
(
Prentice-Hall
,
Englewood Cliffs
,
1971
), p.
194
. Equivalently, one can regard the eigenvalues as the fundamental invariant quantities, corresponding to a maximum of 3 distinct invariants for our 3×3 matrix. Thus, the determinant of Q, (the product of the eigenvalues and an equally attractive prospect for an invariant) can, in fact, be shown to completely determined by the values of TrQ2andTrQ3.
43.

An alternate (and a bit less transparent) derivation is given in Appendix A of Ref. 24.

44.
M. L.
Mehta
, in
Random Matrices
, 2nd ed. (
Academic
,
Boston
,
1991
), Chap. 3.
45.
D.
Chandler
, in
The Liquid State of Matter: Fluids Simple and Complex
, edited by
E. W.
Montroll
and
J. L.
Lebowitz
(
North Holland
,
Amsterdam
,
1982
).
46.
L. E.
Reichl
,
A Modern Course in Statistical Physics
, 4th revised ed. (
Wiley-VCH
,
Weinheim
,
2016
), Appendix A.
47.
Details are provided in  Appendix A.
48.
E.
De Miguel
,
L. F.
Rull
,
M. K.
Chalam
, and
K. E.
Gubbins
,
Mol. Phys.
74
,
405
(
1991
).
49.
L.
Frechette
and
R. M.
Stratt
,
J. Chem. Phys.
144
,
234505
(
2016
).
50.
Details are provided in the supplementary material. See,
Y.
Zhao
, Ph.D. dissertation (
Brown University
,
2018
), for a discussion. The most substantial difference from the protocol published in Ref. 49 is that the gradual compression used to achieve the target densities here starts from a lower density, ρσ3=0.05, and increases the density in smaller, N-dependent, increments ranging from Δρσ3=2×103 (N = 256) to 5×104 (N = 2048).
51.
Matrices were diagonalized via the double-shift Francis method given by
M.
Galassi
,
J.
Davies
,
J.
Theiler
,
B.
Gough
,
G.
Jungman
,
P.
Alken
,
M.
Booth
, and
F.
Rossi
,
GNU Scientific Library Reference Manual
, 3rd ed. (
Network Theory
,
2009
), www.gnu.org/software/gsl. Nematic calculations are carried out in a way that explicitly recognizes that the director direction can fluctuate: the m1,m2, and m3 parameters of Eq. (3.15) are computed with respect to a z axis defined separately for each liquid configuration by that configuration’s principal eigenvector.
52.

It is possible to show from Eq. (3.20) that nematic phases should have ssmacroscopic1/N.

53.
If one defines the correlation length by the requirement that asymptotically 2Cr/Cr=ξ2, then Eq. (4.3) is the (asymptotically zero) function that meets this requirement.
H. E.
Stanley
,
Introduction to Phase Transitions and Critical Phenomena
(
Oxford University
,
Oxford
,
1971
), Chap. 7.
54.
Fits are carried out in the range r/σ = 4–7, with the lower limit needed to avoid the influence of local structure and the upper limit needed to avoid finite-size effects (which are actually visible for r/σ > 7 in Fig. 8). For densities ρσ3 = 0.29, 0.30, 0.304 and the N values used in Fig. 8, the simulation box lengths L/σ = 19.2, 19.0, 18.9, respectively, so the upper limit corresponds to a value between 13L/σ and 12L/σ. Correlation functions computed for smaller N values (not shown) confirm that deviations from the Ornstein-Zernike form at the highest r values are indeed the results of finite-size effects [
Y.
Zhao
, Ph.D. dissertation (
Brown University
,
2018
)].
55.
A.
Widmer-Cooper
,
P.
Harrowell
, and
H.
Fynewever
,
Phys. Rev. Lett.
93
,
135701
(
2004
).
56.
E.
Flenner
and
G.
Szamel
,
Phys. Rev. E
72
,
031508
(
2005
).
57.
M. P.
Allen
and
D.
Frenkel
,
Phys. Rev. Lett.
58
,
1748
(
1987
).
58.

Although single-molecule orientational correlation functions seem to be largely immune to pre-transitional effects in the isotropic phase, they do exhibit a noticeable slowing down once the system crosses into the nematic phase. See Refs. 19 and 23.

59.
L.
Berthier
and
G.
Biroli
,
Rev. Mod. Phys.
83
,
587
(
2011
).
60.
N.
Laevi
,
F. W.
Star
,
T. B.
Schroder
, and
S. C.
Glotzer
,
J. Chem. Phys.
119
,
7372
(
2003
).
61.
C.-P.
Fan
and
M. J.
Stephen
,
Phys. Rev. Lett.
25
,
500
(
1970
).
62.

These divergences are normally described as functions of temperature rather than density, but the Landau analysis is identical when density is the control parameter, as it is for us.

63.
M. D.
Ediger
,
C. A.
Angell
, and
S. R.
Nagel
,
J. Phys. Chem.
100
,
13200
(
1996
).
64.
M. D.
Ediger
and
P.
Harrowell
,
J. Chem. Phys.
137
,
080901
(
2012
).
65.
S.
Ryu
and
R. M.
Stratt
,
J. Phys. Chem. B
108
,
6782
(
2004
).
66.
There are yet other differences between the slowing down of relaxation approaching a glass transition and the slowing down seen during the approach to a weakly first-order transition or critical point. For a nice discussion of this issue and its illustration in the context of hard-disk systems, see,
S.
Yaida
,
L.
Berthier
,
P.
Charbonneau
, and
G.
Tarjus
,
Phys. Rev. E
94
,
032605
(
2016
).

Supplementary Material