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.

## I. INTRODUCTION

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 transition^{2} 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 $\Omega ^j$ (*j* = 1, …, *N*), but which have no net head-or-tail preference, the average value of the second Legendre polynomial $P2x=123x2\u22121$ of the components along the director

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}

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 transition^{29–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 theory^{34} 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 model^{35} 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}

## II. LIQUID CRYSTALLINE ORDER AS A RANDOM MATRIX PROBLEM—THE ISOTROPIC PHASE

### A. Some concepts drawn from random matrix theory

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 spacings^{37} 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,

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

normalized so that

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,

with **O** being an orthogonal matrix).^{41} The natural candidates for such invariants are traces of Q to arbitrary powers, $TrQ,TrQ2,TrQ3,TrQ4,\u2026$, 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, $\rho \lambda 1,\lambda 2,\lambda 3$. These traces themselves are easily expressed in terms of the three eigenvalues $\lambda 1,\lambda 2,\lambda 3$. For example,

Inasmuch as the *Q*_{a} 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 $\lambda 1,\lambda 2,\lambda 3\u2192Q1,\u2026,Q6$, the so-called Vandermonde determinant,^{44}

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* ($=\lambda 1$),^{28,40}

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.

### B. A microscopic statistical mechanical derivation of the extent of orientational order in an isotropic phase

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

tells us that

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 *j* ≠ *k*] 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 *N*_{corr} the typical number of molecules orientationally correlated with a given molecule, then Eqs. (2.8) and (2.9) predict

In typical molecular liquids, we find $Ncorr$ to be on the order of the size of a solvation shell or two^{45} so that $Ncorr/N\u226a1$. 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 $\Omega ^j$ (*j* = 1, …, *N*),

or, using Fourier representations of the delta functions

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*),

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

denote orientational fluctuations of the *j*th molecule.

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

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,

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

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

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

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

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

### C. The implications for finite-size scaling

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 α = (5*N*/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 $\alpha 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

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}

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.

## III. SOME GENERALIZATIONS—NEMATIC PHASES AND TIME DEPENDENCE

### A. Two-time probability distributions

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),\u2009Q(t)$. The list of possible invariants enlarges a bit, but remains short, namely,

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 **Q**^{2} invariants listed in Eq. (3.1).

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

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

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

where, for any *t*,

Using the same Fourier representations of the delta functions,

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

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}

but where, now, $t1\u2009and\u2009t2$ are two, possibly different, times. The generating function therefore continues to take on a simple form

[with $Gk\u2009and\u2009GK$ 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),

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

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,

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

### B. The distribution of the order parameter tensor in the nematic phase

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^\u22c5Q$ and $1\u2212n^n^\u22c5Q$, respectively, it cannot be a function of just the $TrQn$ quantities we had been using. Moreover, since the average $Q$ is no longer zero,

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 $\delta Q=Q\u2212Q$ instead of the **Q** tensors themselves. Thus, the tensor that ends up being involved is the most general tensor that has these properties,

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 $\delta Q2$ invariants in the nematic case,

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

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

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,

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

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

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

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*,

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

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

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),

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 *I*_{3} 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.

## IV. APPLICATION OF RANDOM-MATRIX THEORY TO A LIQUID CRYSTAL SIMULATION

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.

### A. Details of the model and simulation methods

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}

(in the usual notation), and limit ourselves to studying a range of dimensionless densities $\rho \sigma 3$ at a fixed dimensionless temperature,

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 $\delta t=0.001\tau 0$ with $\tau 0=m\sigma 2/\epsilon 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 $\rho \sigma 3=0.304$ were averaged over $1.5\xd7107$ liquid configurations and all of results pertaining to other densities or to uncorrelated systems were averaged over $3\xd7106$ 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

**we sampled. The matrix was then diagonalized by standard methods and the order parameter for that configuration**

*R**s*(

**), 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.**

*R*^{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 $\rho \sigma 3\u22640.304$ correspond to single-phase macroscopically isotropic liquids, while those with reduced densities $\rho \sigma 3\u22650.33$ correspond to single-phase macroscopically nematic liquids.

### B. Random-matrix theory and finite-size scaling

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.

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).

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 $\rho \sigma 3\u22640.30$, and densities $\rho \sigma 3>0.30$ are precisely where we see the random matrix formulas begin to converge slowly.

### C. Random matrices as a window into pseudonematic domains

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

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

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}

. | Single-molecule properties^{b}
. | Collective properties^{c}
. | |||
---|---|---|---|---|---|

ρσ^{3}
. | $\gamma 1st-shell$ . | $\gamma $ . | $\tau /\tau 0$ . | $\xi /\sigma $ . | $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 properties^{b}
. | Collective properties^{c}
. | |||
---|---|---|---|---|---|

ρσ^{3}
. | $\gamma 1st-shell$ . | $\gamma $ . | $\tau /\tau 0$ . | $\xi /\sigma $ . | $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/\epsilon $ = 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\xd7107$ configurations for susceptibilities, and over $1.2\xd7107$ and $5\xd7104$ 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 $\tau 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=\rho \xi 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,

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

Hence, the dimensionless non-Gaussian parameter

should vanish as the system size increases.

As is evident from Fig. 9, before the intermediate density regime is reached, this $\alpha 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 $\rho \sigma 3=0.304$, a density in the intermediate regime, $\alpha 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.

We suggest that the key to understanding this dichotomy is to realize that this threshold for Gaussian behavior is of the order of $Ncorr=\rho \xi 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),

will always cause the two-particle invariant to dominate, leading to a Gaussian distribution of orientational order. However, if $N\u2264Ncorr$, 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).

### D. Dynamical implications of pseudonematic domains

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),

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}

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=\gamma t/\gamma $, 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 $\rho \sigma 3=0.304$.^{58}

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,

to compute the corresponding random-matrix-theory (RMT) averages at the relevant RMT values: $A=\Gamma (t),B=C(t)\Gamma (t)$. Since $ZA,B\u223c\alpha \u22125A2\u2212B2\u22125/2$,

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,

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 $\alpha 2t$,

Since *C*(0) = 1, it is easy to confirm that $\alpha 20=\alpha 2$, the static non-Gaussian parameter of Eq. (4.4b), and since $\alpha 2\u221e=0$, the decay of $\alpha 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 $\alpha 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

precisely the kind of correlations crucial to the $\chi 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 $\rho c$,

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),

which itself is simply related to the correlation length. Integrating the asymptotic form, Eq. (4.3), for example, would imply that $\gamma =4\pi \xi 2$ and hence that we should expect

## V. CONCLUDING REMARKS

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 crystals^{10–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 *P*_{2} 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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.

### APPENDIX A: ORIENTATIONAL AVERAGES IN THE ISOTROPIC PHASE

Both the developments in Sec. II and in Sec. III A require that we evaluate averages over the single-molecule orientational tensor elements $q\mu \nu j$ [introduced in Eq. (2.13)] in isotropic systems $(j=1,\u2026,N;\mu ,\nu =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 $\Omega ^j$ is a unit vector, $\u2211\mu \Omega j\mu 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 $\Omega j\mu \u2192\u2212\Omega j\mu $, in an isotropic system, we have

Hence, all single-molecule averages must vanish,

Similarly, averages over the product of two tensor elements $q\mu \nu jt1q\gamma \delta kt2$ (involving what may be two different molecules *j* and *k* at two different times $t1andt2$) can be evaluated using the isotropic-system identities

Equation (A3) follows from the unit-vector property $\u2211\mu ,\nu \Omega j\mu 2t1\Omega k\nu 2t2=1$, and Eq. (A4) can be derived by noting that averages must be invariant to 45° rotations around any third Cartesian axis $\gamma \u2260\mu ,\nu $,

for any *i* and *t*.

Inasmuch as we can take an axis *μ* to point in any arbitrary direction $\Omega ^$,

and that all other averages of the form $q\mu \nu jt1q\gamma \delta kt2$ vanish.

### APPENDIX B: ORIENTATIONAL AVERAGES AND INVARIANTS IN THE NEMATIC PHASE

The expressions derived in Appendix A for averages over the single-molecule orientational tensor elements $q\mu \nu 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 $\Omega j\mu \u2192\u2212\Omega j\mu $. Suppose we consider the consequences of these residual symmetries, calling *μ* = *z* the director axis.

As before, the components of each orientational vector $\Omega ^j$ satisfy the normalization condition, $\Omega jx2+\Omega jy2+\Omega 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

so that

All other single-molecule averages must vanish

Averages over pairs of tensor fluctuations

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 $Tr\u2009Q2$ 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 $\delta qaj\delta 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 *x* ↔ *y*) 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,\u2009m2,\u2009m3,\u2009m4$ defined in Eq. (3.15). In terms of the particular coordinate system we are using,

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,

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

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 *dφ* about the director $n^$, the orientation of the *j*th molecule transforms $\Omega ^j\u2192\Omega ^j+d\phi n^\xd7\Omega ^j$ so that fluctuations in the orientational order parameter tensor,

themselves transform as

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^$,

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

any two-body invariant that has components taken from $\delta 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,

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

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

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 *x* ↔ *y* invariance and Eq. (B4), we obtain

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

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

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,

and thus, in the isotropic limit,

## REFERENCES

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

*p*th order polynomial equation), so any $Qn$ (

*n*> 3) can be expressed as a linear combination of $Q,Q2,and\u2009Q3$. See, for example,

**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$.

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

*N*-dependent, increments ranging from $\Delta \rho \sigma 3=2\xd710\u22123$ (

*N*= 256) to $5\xd710\u22124$ (

*N*= 2048).

*z*axis defined separately for each liquid configuration by that configuration’s principal eigenvector.

It is possible to show from Eq. (3.20) that nematic phases should have $s\u2212smacroscopic\u223c1/N$.

*N*values used in Fig. 8, the simulation box lengths $L/\sigma $ = 19.2, 19.0, 18.9, respectively, so the upper limit corresponds to a value between $13L/\sigma $ and $12L/\sigma $. 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 [

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.