# The nematic-isotropic transition of the Lebwohl–Lasher model revisited

## Abstract

We revisited the nematic-isotropic (NI) transition of the Lebwohl–Lasher lattice model with a detailed investigation of samples containing 200 × 200 × 200 particles. The large-scale Monte Carlo (MC) simulations involved were carried out following the standard Metropolis, as well as the cluster MC Wolff algorithms. A notable free-energy barrier was observed between the isotropic and nematic phase, leading to long-lived metastable states and hysteresis. We provide an improved estimate of the nematic-isotropic transition temperature *T*_{NI}, of the supercooling and superheating temperatures, of the temperature of divergence of pretransitional effects ${T}_{C}^{\ast}$ as well as an analysis of the size distribution of the ordered domains above *T*_{NI}, contributing to a better understanding of this transition of key importance for liquid crystals.

This article is part of the theme issue ‘Topics in mathematical design of complex materials’.

### 1. Introduction

Lattice models, where interaction centres (spins) with *n*_{o} degrees of freedom occupy the sites of a certain *n*_{d}-dimensional lattice, have played and play an extremely important role in condensed matter physics, albeit being apparently removed from the physical systems (magnets, liquids, liquid crystals) they try to address. Indeed, notwithstanding their simplicity, successful lattice models capture the essential symmetries of the problem they apply to, allowing, in the absence of analytic solutions (except of course for the immensely important *n*_{o} = 1, *n*_{d} = 2 Ising model), the in-depth study of issues like phase transitions and defects. The role of lattice models is also, in connection with numerical approaches like Metropolis Monte Carlo (MC) [1], that of assisting the development of methodologies, like cluster moves or histogram reweighting [2,3] in MC that can then be employed for other much more complex systems. In the field of liquid crystals (LC) [4], one of the fundamental issues still debated is the nature of the isotropic-nematic phase transition, occurring at a certain temperature *T*_{NI}, which is experimentally observed to be a weak first-order one with very important pretransitional effects observed in light scattering [5], and in electric field (*Kerr effect*) [6] or magnetic field (*Cotton-Mouton effect*) [7] induced ordering measured, e.g.by the induced birefringence Δ*n* above *T*_{NI}. These effects are reminiscent of a second-order transition and indeed they are supposed to be manifestations of a virtual second-order transition occurring at a temperature *T**_{C} slightly below *T*_{NI} [8,9]. A rationalization of pretransitional effects based on the Landau-de Gennes theory [4] predicts a universal behaviour for nematogens, with exponent $\beta =\frac{1}{2}$ for the orientational order parameter 〈*P*_{2}〉. Early experiments on *n-(4-methoxybenzylidene)-4-butylaniline* (MBBA) seemed to comply with this, except for a critical region very close to the phase transition [5]. Similar behaviour was observed also for *4′- n-pentyl-4-cyano biphenyl* (5CB) from Kerr Effect [10] and NMR [7] measurements. However, other experiments and analyses have suggested instead a tricritical behaviour and exponent

*β*= 0.245 [11]. More recent works both on MBBA and other families of nematogens [12] have found different results for

*T**

_{C}−

*T*

_{NI}and non-universal exponents, also criticizing the lack of adequate chemical stability of previously studied nematogens like MBBA. One possible perturbing effect in experiments might be the existence of a smectic transition in some of the materials studied (e.g.

*4′-*(8CB)) and, in general, the importance of chemical details such as molecular biaxiality and/or flexibility. To eliminate these spurious effects, and to see to what extent the essential features of the NI transition, like pretransitional effects, the existence of a

*n*-octyl-4-cyano biphenyl*T**

_{C}and its closeness to

*T*

_{NI}, are maintained, irrespective of details, the computer simulations of simple lattice models where constituent particles (spins) are truly uniaxial and no smectic or crystal phase exist, appear to be essential.

The reference lattice model in the LC field is the Lebwohl–Lasher (LL) one [13] that was proposed nearly 50 years ago to study the orientational transition between the nematic and isotropic phase. Similarly to the famous *n*_{o} = 3, *n*_{d} = 3 ferromagnetic Heisenberg model [14–17], put forward to study magnetic effects in solids, it considers an array of uniaxial spins placed at the sites of a lattice (simple cubic in the LL case) and free to rotate in 3D, while interacting only with their nearest neighbours. In the Heisenberg model, the demagnetization transition is second-order and the reduced critical temperature, estimated with extensive MC work [15], is *T** = *k*_{B}*T*/*ϵ* = 1.442929(4) (*T* here denotes the absolute temperature, *k*_{B} is the Boltzmann constant, and *ϵ* the characteristic energy scale), while the critical exponents correspond to the behaviour of isotropic ferromagnets such as Fe, Ni and EuO [18]. Unlike the polar spins of the Heisenberg model, the LL spins, intended to model molecules (or tightly packed clusters of molecules [19]) forming non-polar nematics, are non-polar. More specifically, the LL interaction energy for a pair of nearest-neighbour particles whose long axes are represented by headless unit vectors **u**_{i} and **u**_{j} is given by

*ϵ*

_{ij}is a characteristic interaction energy equal to a positive constant,

*ϵ*, when

*i, j*are nearest neigbours and 0 otherwise, while

*P*

_{2}(

*x*) ≡ (3

*x*

^{2}− 1)/2 is the second Legendre polynomial. The cubic lattice spacing,

*a*, is related to molecular size and is typically of the order of a few nm.

The LL model has a weakly first-order orientational transition, like the nematic-isotropic (NI) one in real LC systems. Thus the situation is more difficult to analyse than that of continuous transitions, as we do not have the tools of universality class and critical exponents. However, a helpful analysis has been proposed and applied to *q*-state Potts models that change from second to first-order transition as the integer *q* labelling the number of allowed spin orientations increases to *q* = 5 and above [20–23].

As there is no general analytical solution of the LL model—and, differently from Ising or Heisenberg models, only the first couple of non-trivial terms in the low- and high-temperature series expansion of the free energy have been evaluated [13]—computer simulations have been used to investigate the nematic-isotropic transition. Already the early Metropolis MC [1,24] simulations have shown [13], as already mentioned, that there is a weakly first-order phase transition [25] between the orientationally ordered nematic phase, where there is an average spin orientation given by a unit vector called the director, and the disordered isotropic phase. The early LL work was followed by further Metropolis MC studies that gradually increased the system size and improved the NI transition temperature estimate [26,27] and later also by advanced simulation techniques reaching beyond Metropolis [3,28–30].

Inevitably, all of these studies have dealt with a finite spin system, which strictly suppresses the first-order (discontinuous) character of the NI transition, even if periodic [31] or cluster [27] boundary conditions are applied at the simulation box boundaries to mimic a larger bulk sample. Finite-size analyses performed by Zhang *et al.* [28] and later by Priezjev & Pelcovits [29] have shown that the first-order character of the transition is gradually restored when the system size is increased (up to 70^{3} spins in the latter case), and provided an estimate of the transition temperature for an infinite sample (i.e. in the thermodynamic limit). However, it is worth noting that finite-size scaling methods [25] are of less straightforward application for temperature-driven first order transitions as, for instance, shown by Binder and others [21,23,32] in their study of the *n*_{d} = 2, *q*-states Potts models where lattices as large as 240 × 240 have been deemed necessary [23].

In the present paper, almost half a century after the pioneering work by Lebwohl & Lasher [13], and after a significant increase in computing power, we consider an even larger sample containing *N* = 200^{3} spins with periodic boundary conditions (PBC), whose linear size well exceeds the nematic-isotropic correlation length (estimated to be of the order of ∼10 *a* at the NI transition) by more than one order of magnitude and thus, we expect to observe genuine first-order NI behaviour in MC simulations. The aim of the work is to exploit the results to tackle a number of open issues related to the NI transition.

The first is refining the accuracy in the estimate of the transition temperature. The most recent values for *T*_{NI} obtained, with various methods, from [33] and [30] differ by up to 4 × 10^{−4} in reduced temperature scale, ${T}^{\ast}={k}_{B}T/\u03f5$, which, scaled to the transition temperature of an ordinary room-temperature nematic LC, would correspond to more than 0.1 K, a very large uncertainty considering that real experiments aimed at studying transition details are some two orders of magnitude more precise, at mK level [34].

In addition, we would like to refine the estimate for the divergence temperature of pretransitional effects, ${T}_{C}^{\ast}$, that 25 years ago in the 30 × 30 × 30 MC simulations of [26] was found to be about 1 K equivalent below *T*_{NI} in agreement with the experiment, and see if that estimate still holds for a sample two orders of magnitude bigger.

Another point of interest is to see if there is a change of pretransitional behaviour on approaching the transition from above, as the temperature is lowered towards the weakly first-order transition temperature [5,23].

In the next section, we describe the simulation performed and discuss the results and data analysis in detail, following which we summarize our main findings.

### 2. Monte Carlo simulations and results

In a first attempt, we performed standard Metropolis MC simulations involving individual spin moves. Every MC cycle consisted of an attempt to reorient each of the spins individually, following a spin sequence obtained from a random shuffling algorithm. The trial spin orientations were generated following the Barker–Watts technique [35] and were automatically accepted if the energy of the new configuration (*U*_{n}) was lower than the energy of the old one (*U*_{o}), or accepted with probability exp[*β*_{T}(*U*_{o} − *U*_{n})] otherwise [1]; here *β*_{T} ≡ 1/(*k*_{B}*T*). The trial reorientation amplitude was adjusted dynamically in each cycle to maintain an acceptance ratio of around 50%. After equilibration, several averages of interest were calculated.

The reduced interaction energy per lattice site is given by ${U}^{\ast}={(N\u03f5)}^{-1}\sum _{\u27e8i,j\u27e9}{U}_{ij}$ and its average over 1.2 × 10^{6} MC production cycles (1 cycle corresponding to *N* attempted individual MC moves) is plotted versus reduced temperature ${T}^{\ast}\equiv {k}_{B}T/\u03f5$ in figure 1*a*, with a close-up of the transition region in figure 1*b*. The plotted results were obtained by launching independent Metropolis MC runs at selected temperatures either from isotropic or perfectly aligned initial conditions, the experimental equivalents of this approach being cooling and heating temperature scans, respectively. These datasets are shown together with a third dataset obtained from simulations based on the Wolff cluster MC technique implementing collective rather than individual spin reorientations [36], started from the isotropic phase in our case. We can see excellent agreement from all three sets of results and that cooling down and heating up results are perfectly superimposable, confirming proper equilibration, except in the very narrow temperature range close to the NI transition, where hysteresis, i.e. supercooling and respectively superheating from runs started from isotropic or ordered states appears (as apparent in the close-up shown in figure 1*b*). Hysteresis is typical of first-order transitions and experimentally observed. For our very large samples (and differently from the sample of a few tens of thousand particles [26]), it was obtained systematically when launching the simulation from aligned nematic and disordered isotropic states, respectively.

The orientational order was monitored by diagonalizing the second-rank ordering matrix $\mathbf{\text{Q}}={N}^{-1}\sum _{i=1}^{N}(3{\mathbf{\text{u}}}_{i}\otimes {\mathbf{\text{u}}}_{i}-\mathsf{I})/2$ in a given MC cycle; then, the largest eigenvalue of **Q** was averaged appropriately over MC cycles to give the orientational order parameter 〈*P*_{2}〉 shown in figure 2*a*,*b*. In the isotropic phase one has 〈*P*_{2}〉 ∼ 0, while a non-zero 〈*P*_{2}〉 is a signature of nematic order. Note that the above algorithm always gives a positive contribution to the 〈*P*_{2}〉 average and therefore the result obtained in the isotropic phase inevitably differs from zero (even if not significantly). The middle eigenvalue of **Q** denoted by $\u27e8{P}_{2}^{m}\u27e9$, also shown in figure 2*a*,*b*, does not suffer from this problem [37] and can be used to more reliably detect the absence of orientational ordering in the isotropic phase.

The order parameters 〈*P*_{2}〉, $\u27e8{P}_{2}^{m}\u27e9$ show a discontinuous jump at the NI transition. A more detailed look at the transition region (figure 2*b*) reveals, like for the energy in figure 1*b*, the existence of hysteresis and long-lived metastable states—superheated nematic and supercooled isotropic phase.

The heat capacity per site (and in units of *k*_{B}) can be obtained from interaction energy fluctuations as ${c}_{V}^{\ast}=N[\u27e8{({U}^{\ast})}^{2}\u27e9-{\u27e8{U}^{\ast}\u27e9}^{2}]/{({T}^{\ast})}^{2}$. Alternatively, it can also be calculated directly from its definition, ${c}_{V}^{\ast}=\text{d}\u27e8{U}^{\ast}\u27e9/\text{d}{T}^{\ast}$, yielding similar results. In our analysis, the latter approach was used, approximating the derivative by a finite difference, and the results, shown in figure 3*a*,*b*, exhibit a pronounced ${c}_{V}^{\ast}$ peak in the transition region.

In the temperature range between ${T}^{\ast}\sim 1.12220$ and ${T}^{\ast}\sim 1.12255,$ we simultaneously observe persistent nematic and isotropic states at a given temperature, and no switching between them during the 1.2 × 10^{6} MC production cycles is observed. While for the relatively small system sizes of *N* = 30^{3} at the NI transition the system frequently switches between the two states and the transition temperature can be determined from e.g. order parameter histograms by comparing the nematic and isotropic peak magnitudes [26]; in the current much larger system this is made impossible by the much higher free energy barrier separating the nematic and isotropic states, known to increase with system size [38]. This problem has been overcome in [29] by performing collective, rather than individual, spin reorientations, following the Wolff cluster algorithm [36]. While this solution works well for *N* = 70^{3} considered in [29] and is to some extent effective even for some *N* = 100^{3} runs we performed here, no hopping between the nematic and isotropic free energy minima could be observed during the production run for the current system size, *N* = 200^{3}. The absence of hopping between the disordered and ordered basins is compatible with the observed energy and order parameter histograms where the nematic and isotropic branches are never even close to overlapping, which makes a determination of NI transition temperature *T*_{NI} in this manner impossible.

An alternative route to determine *T*_{NI} would be to follow the original LL work [13], making use of thermodynamic integration. As already mentioned, in the low- and high-temperature limits, explicit expressions for the free energy ${F}^{\ast}({\beta}_{T})$ can be obtained. Moreover, at any temperature one can obtain the average energy $\u27e8{U}^{\ast}\u27e9$ from simulation, and integrate the relation $\u27e8{U}^{\ast}\u27e9=\text{d}({\beta}_{T}{F}^{\ast})/\text{d}{\beta}_{T}$ over *β*_{T} to reconstruct the temperature dependence of ${F}^{\ast}$ for both branches, nematic and isotropic,

*T*

_{NI}value. However, it turns out that the required four-digit precision cannot be achieved in this way. The procedure is quite sensitive to the choice of the high- and the low-temperature reference, number and position of data points used in the integration, the accuracy of the corresponding average energy $\u27e8{U}^{\ast}\u27e9({\beta}_{T})$ representing the slope of ${\beta}_{T}{F}^{\ast}$ when plotted versus

*β*

_{T}, and even on the numerical integration method.

Another possibility to find *T*_{NI} is to attempt a ‘coexistence’ simulated experiment wherein the starting configuration for an MC simulation at a given temperature is assembled by taking, say, the bottom half of the sample from the nematic branch and the upper half from the isotropic branch, both pre-equilibrated at that temperature. The Metropolis MC evolution (simulating the canonical ensemble and thus minimizing the total free energy) should then eventually lead to the equilibrium, either nematic or isotropic, at the given temperature. We have performed a number of runs involving different spatial partitions of the nematic and isotropic phase in the initial sample (e.g. layered or checkerboard-like in 3D, with different domain sizes), as well as employing different random-number generators, starting configurations, and random seeds, using either the Metropolis or Wolff algorithm. In these experiments, it was essential that the initial nematic and isotropic volumes be exactly the same due to the high surface-to-volume ratio of the thus generated domains to avoid bias and possible issues with interfacial tension. The random-number generators used were the xoshiro256** [39], a recent high-quality generator using combinations of linear transformations and scramblers (chosen to be the default one in our study), and the standard ran2 generator from Numerical Recipes [40] (for a quick check and comparison). After performing a total of 24 runs at each temperature value, all runs for ${T}^{\ast}\le 1.12235$ turned nematic, while all runs for ${T}^{\ast}\ge 1.12250$ turned isotropic. The runs at intermediate values (1.12240 and 1.12245) turned either nematic or isotropic, subject to thermal noise. From these observations, we estimate the NI transition temperature in an *N* = 200^{3} sample, within four-digit precision, as ${T}_{NI}^{\ast}={1.1224}_{-0.00003}^{+0.00007}$. Moreover, from figures 1*b*–2*b*, we deduce that the isotropic phase can be supercooled at least down to ${T}^{\ast}=1.12220$, while the nematic can be superheated at least up to ${T}^{\ast}=1.12255$. The actual metastability range is probably somewhat broader since, in the ‘coexistence’ experiment performed close to the two limiting temperatures, the system could hop across the relatively low barrier and escape into the other (deeper) free-energy basin.

Now the transition enthalpy (i.e. the latent heat of the NI transition) Δ*H* can be obtained by taking the internal energy difference of the two coexisting phases at ${T}_{NI}^{\ast}$, yielding $\mathrm{\Delta}{H}^{\ast}=\mathrm{\Delta}H/\u03f5\sim 0.065$. (Reference [28] here gives the very different value Δ*H*/*ϵ* ∼ 0.20 ± 0.04 from finite-size scaling behaviour considering samples up to 28 × 28 × 28 spins.) Similarly, the transition entropy can be estimated as $\mathrm{\Delta}{S}^{\ast}=\mathrm{\Delta}{H}^{\ast}/{T}_{NI}^{\ast}\sim 0.058$. ([26] gives $\mathrm{\Delta}{S}^{\ast}\le 0.05$ considering 30 × 30 × 30 spins.)

Since the latent heat at the weakly first-order NI transition is rather small, various physical observables might be expected to exhibit critical power-law behaviour similar to that observed in second-order transitions, at least very close to the transition [41]. We have thus attempted to estimate the corresponding critical exponents *β, α* and *ν* [42] by fitting the temperature dependencies of the order parameter, specific heat, and correlation length, respectively. The fits turn out not to be very robust and are dependent on the fitting range as in real experiments [41]. On the one hand, the fitting temperature range should be narrow to remain in the vicinity of the NI transition and, at the same time, broad enough to exploit as many data points as possible for a representative fit. Considering the order parameter 〈*P*_{2}〉 in the nematic phase, see figure 2*a*,*b*, its temperature dependence can be fitted to the four-parameter functional form $A+B{({{T}_{C}^{\ast}}^{\prime}-{T}^{\ast})}^{\beta}$. Taking all the data points between ${T}^{\ast}=1.1000$ and 1.12255 (corresponding to ∼7 K for a room-temperature nematic, which is more than in a typical experiment) yields *A* ∼ 0.13, *B* ∼ 0.9, ${{T}_{C}^{\ast}}^{\prime}\sim 1.1226$ and *β* ∼ 0.3. (Experiments give here *β* ∼ 0.25 [41].) Similarly, the temperature dependence of the heat capacity ${c}_{V}^{\ast}$, figure 3*a*,*b*, can be fitted to the three-parameter form ${C}^{\prime}{({{T}_{C}^{\ast}}^{\prime}-{T}^{\ast})}^{-{\alpha}^{\prime}}$, yielding (in the same temperature interval as for 〈*P*_{2}〉), *C*′ ∼ 1.3, ${{T}_{C}^{\ast}}^{\prime}\sim 1.1226$, and *α*′ ∼ 0.4, the experimental values for *α*′ usually falling between 0.3 and 0.4 [41]. Above the NI transition, the form $C{({T}^{\ast}-{T}_{C}^{\ast})}^{-\alpha}$ fits the temperature dependence of ${c}_{V}^{\ast}$ across the temperature interval between ${T}^{\ast}=1.1222$ and 1.1500 (corresponding to ∼8 K) for *C* ∼ 0.5, ${T}_{C}^{\ast}\sim 1.1219$ and *α* ∼ 0.4, while an exponent between 0.11 and 0.5 is expected here [41].

The above results can be compared with the widely used Landau theory [4] based on the phenomenological expansion of free energy *F* into powers of the order parameter, 〈*P*_{2}〉, up to fourth order:

Here *a, b*, and *c* are positive constants, ${T}_{\ast}$ is the lowest supercooling temperature of the isotropic phase, and *F*_{0}(*T*) is a reference free energy that depends only on temperature. In this framework, the NI transition takes place at ${T}_{NI}={T}_{\ast}+2{b}^{2}/9ac$, while the superheating temperature limit of the nematic phase is reached at ${T}^{\u2020}={T}_{\ast}+{b}^{2}/4ac$ [4]. Furthermore, the theory predicts a critical exponent *β* = 1/2 for the order parameter and *α*′ = 1/2, *α* = 0 for the specific heat below and above the NI transition, respectively. While our estimates of these exponents are not extremely reliable, however, they do clearly differ from the predictions of the Landau theory. Moreover, in the Landau theory, the ratio of the temperature differences between the NI transition and the supercooling limit of the isotropic phase and between the superheating limit of the nematic and the NI transition equals $({T}_{NI}-{T}_{\ast})/({T}^{\u2020}-{T}_{NI})=8$, while from our simulated results it is roughly estimated as ∼1.33, which is very different. Alternatively, the fitted parameters ${T}_{C}^{\ast}$ and ${{T}_{C}^{\ast}}^{\prime}$ give a ratio of ∼2.5, also different from the Landau theory prediction. All these differences between Landau theory and MC predictions of the LL model should not be surprising—the former is a mean-field approach entirely neglecting fluctuations even close to the transition, while the latter are based on pairwise interactions and offer a more adequate description of critical phenomena.

It is also interesting to study orientational correlations in the LL model system. In particular, the correlation length *ξ*, which remains finite at the first-order NI transition, should be compared with the linear size of the simulation box to assess finite-size effects. To analyse correlations between spins separated by a distance *r*, the orientational correlation function can be conveniently introduced as *G*_{2}(*r*) = 〈*P*_{2}(**u**_{i} · **u**_{j})〉_{r} [26], where the average 〈 · · · 〉_{r} runs over all eligible spin pairs (*i, j*), and over MC cycles. To reduce the computational effort, every 10 MC cycles we selected 100 spins at random and used them in the calculation of *G*_{2}(*r*). At each temperature, the resulting *G*_{2}(*r*) were fitted to the Ornstein–Zernike form *G*_{2}(*r*) = (*D*/*r*)exp( − *r*/*ξ*), up to a cut-off at *r* ∼ 85 *a*; here *D* and the correlation length *ξ* are the fitting parameters [26]. The temperature dependence of *ξ* in the isotropic phase is shown in figure 4; it is largest at the supercooling limit of the isotropic phase where it reaches ∼18 *a*, and is equal to ∼15 *a* at the NI transition. This is more than one order of magnitude below the linear dimension of the simulated sample, 200 *a*, therefore, we presume that the features of the observed (weakly) first-order NI transition are already close to those in the thermodynamic limit, including the transition temperature. The temperature dependence of *ξ* itself can also be fitted to the functional form $E{({T}^{\ast}-{T}_{C}^{\ast})}^{-\nu}$, yielding *E* ∼ 0.6, ${T}_{C}^{\ast}\sim 1.1219$, and a critical exponent *ν* ∼ 0.4 when taking all the data points shown in figure 4. Note that here the Landau theory prediction is *ν* = 1/2.

An alternative way of determining the divergence temperature of pretransitional effects, ${T}_{C}^{\ast}$, is to examine the orientational correlations in the isotropic phase by calculating the so-called second-rank Kirkwood coefficient *g*_{2} [26,43] that diverges at a second-order transition. For lattice simulation purposes, it is defined by ${g}_{2}=\sum _{k}{z}_{k}{G}_{2}({r}_{k})$, where *k* runs over spherical shells of radius *r*_{k}, each of them containing *z*_{k} *k*-th nearest neighbours of a selected particle on the lattice. Up to half the simulation box, *r*_{k} = 100 *a, g*_{2} can be calculated directly from the simulated *G*_{2}(*r*) for each temperature, but close to the NI transition the sum over *k* still does not fully converge. Simply extending the sum to larger *r*_{k} is impossible due to the PBC applied; instead, the simulated *G*_{2}(*r*) can be replaced by their Ornstein–Zernike fits used to sufficently extend the summation range [26], in our case up to a cut-off at *r*_{k} = 400 *a*. The Landau theory now predicts a linear dependence when plotting ${T}^{\ast}/{g}_{2}$ versus ${T}^{\ast}$, with the intercept giving the value of ${T}_{C}^{\ast}$ [43]. The corresponding dependence plotted with our data is shown in figure 5*a*,*b* and it turns out that close to the NI transition, there is a deviation from the predicted linear behaviour, similar to what is observed experimentally [12]. When taking data points for $1.1260\le {T}^{\ast}\le 1.1500$, i.e. not too close to the transition where the linear dependence is maintained over a relatively large temperature range and the Landau theory appears to remain valid, linear extrapolation (line A) yields an intercept at ${({T}_{C}^{\ast})}_{A}=1.1202$, which is in excellent agreement with the result obtained in [26]. On the other hand, taking the data close to the NI transition in the rather narrow range $1.1222\le {T}^{\ast}\le 1.1225$ results in an intercept at ${({T}_{C}^{\ast})}_{B}=1.1219$ (line B), which, on the other hand, agrees with the estimates obtained from the fits of ${c}_{V}^{\ast}$ and ${\xi}^{\ast}$. Moreover, this ${T}_{C}^{\ast}$ value results in ${T}_{NI}^{\ast}-{T}_{C}^{\ast}\sim 5\times {10}^{-4}$, corresponding to ∼0.15 K in real units, while the extrapolation through line A would give ∼0.7 K, which is quite different (figure 5*b*).

Motivated by the two different fits, we tried again to fit the temperature dependence of the correlation length ${\xi}^{\ast}({T}^{\ast})$, considering only those temperatures not too close to ${T}_{NI}^{\ast}$ that resulted in line A in figure 5, and saw a change of the critical exponent *ν* from ∼0.4 reported earlier to ∼0.5. This again agrees well with the simulation result of [26]. Note that the value *ν* ∼ 0.5 agrees also with the prediction of the Landau theory, which, together with the limited linearity in the ${T}^{\ast}/{g}_{2}$ versus ${T}^{\ast}$ plot suggests that the Landau description may be suitable in this temperature range, but not closer to the transition.

In pretransitional ordering, the orientational domains can also be visualized directly. In a selected snapshot, for each lattice site the ordering matrix **Q**_{i} was averaged over lattice sites up to a distance of $a\sqrt{5}$, which gave a total of 57 spins per site. After diagonalizing **Q**_{i}, the eigenvector belonging to the largest eigenvalue was identified as the local nematic director **n**_{i} and colour-coded according to its orientation, while its contrast was taken proportional to the degree of order quantified by the first Westin metric [44] element, *c*_{l} in the notation of Callan-Jones *et al*. [45]. (A similar method has been applied to visualize the orientational domains in polydomain liquid crystal elastomers [46].) Examples of resulting coarse-grained snapshots in the isotropic phase are shown in figure 6. Cases (*a*) and (*b*) correspond to the supercooling limit at ${T}^{\ast}=1.1222$ and to ${T}^{\ast}=1.1500$ deep in the isotropic phase, respectively. The difference in orientational domain sizes between the two cases is visible with the bare eye — the domains are larger at the supercooling limit. Next to each snapshot, sample cross-sections are shown to depict the domains oriented approximately in the same direction, e.g. along the *z*-axis. While deep in the isotropic phase the domains appear to be mostly isolated droplets, at the supercooling limit they become larger and take more complex shapes, extending quite far across the sample, which is compatible with the previously observed increase in correlation length *ξ*.

For a more quantitative analysis of droplet size and shape, we calculated the distribution of domain sizes, following the approach of [46]. In the procedure, an orientational domain is assigned to each lattice site *i*. The domain size is determined from a recursive walk on the cubic lattice starting from site *i* and identifying all lattice sites *j* in its vicinity with directors **n**_{j} that satisfy *P*_{2}(**n**_{i} · **n**_{j}) > 0.85 and *c*_{l} > 0.12—i.e. counting well-ordered sites with the local director almost parallel to the initial **n**_{i}. Note that, unlike in models with discrete states like Ising or Potts, here the criterion for domain identification is somewhat arbitrary. The domains defined in the above manner heavily overlap, however, looking at the statistics of their sizes (i.e. the number of sites they contain) still provides some insight into the spatio-orientational organization of the system. Figure 7 shows the normalized domain size histogram *p*(*s*) plotted against the domain size *s*. At the supercooling limit, as well as not too far from the NI transition, the log-log plot reveals a linear dependence compatible with a power law $p(s)\propto {s}^{\tau}$, where *τ* ∼ −1.4, over a broad (but not entire) range of domain sizes *s*. This corresponds to, e.g. case (*a*) in figure 6 and may be a signature of fractal domains in the system [23] characteristic of a second-order transition. On the other hand, away from the transition, case (*b*) in figure 6, strong deviations from linear behaviour in the log-log plot are observed. Strictly speaking, at the transition itself (dashed curve in figure 7), the linear dependence breaks down at *s* ∼ 10^{4}, while at the supercooling limit it persists up to a somewhat larger value. The *p*(*s*) dependence at *T*_{NI} is qualitatively similar to that obtained for the *q* = 5 Potts model [23] from where we can again conclude that the NI transition in the LL model is first-order, however, with more and more pronounced features of a second-order transition (i.e. extended linear range in figure 7 and an additional increase of *ξ*) when cooling down in the metastable isotropic branch towards the supercooling limit and ${T}_{C}^{\ast}$.

The transition temperature observed in a finite system has been predicted to decrease with increasing system size [38], which has also been seen in several LL simulation studies [28–30]. There, authors have performed finite-size scaling analyses extrapolating their ${T}_{NI}^{\ast}$ values obtained in smaller samples to *N* → ∞, the thermodynamic limit. To obtain ${T}_{NI}^{\ast}$ for a given system size, they typically analysed the relative depth of the nematic and isotropic free-energy minima and/or the peak behaviour of the specific heat and orientational susceptibility [28–30]. Zhang *et al.* [28] studied systems of up to *N* = 28^{3} and predicted ${T}_{NI}^{\ast}=1.1232\pm 0.0001$, which agrees with the early analysis by Fabbri and Zannoni in an *N* = 30^{3} sample [26]. Furthermore, Shekhar *et al.* [30] considered lattices up to *N* = 40^{3} and give ${T}_{NI}^{\ast}=1.1226\pm 0.0001$, 1.1228 ± 0.00014 or 1.1229 ± 0.00015, depending on the ${T}_{NI}^{\ast}$ determination method. Largest samples so far, up to *N* = 70^{3}, have been analysed by Priezjev & Pelcovits [29] who obtained ${T}_{NI}^{\ast}=1.1225\pm 0.0001$ for *N* → ∞. Results from [26,29,30] are summarized in figure 8, together with our current estimate ${T}_{NI}^{\ast}=1.1224{\hspace{0.17em}}_{-0.00003}^{+0.00007}$ for the transition temperature in a finite but large (*N* = 200^{3}) sample. Our estimate lies below all cited ${T}_{NI}^{\ast}$ values but is compatible (within error bars) with the *N* → ∞ prediction by Priezjev and Pelcovits, while Shekhar *et al*. appear to overestimate the transition temperature regardless of the method employed.

### 3. Conclusion

We have revisited the nematic-isotropic transition of the Lebwohl–Lasher model, performing a large-scale Monte Carlo simulation in a model cubic lattice system with periodic boundary conditions whose linear size (200 lattice spacings) exceeds the orientational correlation length observed at the transition (approx. 15 lattice spacings) by more than one order of magnitude. Dealing with these system sizes has highlighted the difficulties, related to persistent metastable states and hysteresis of this weak first-order transition, in locating the NI transition with techniques such as energy and order parameters histograms analysis [26] and reweighting [29,30] successfully employed for samples one to two orders of magnitude smaller. After testing various additional approaches, as described in the previous section, we have refined the estimate for the transition temperature, along with an accurate determination of its superheating and supercooling bounds, to ${T}_{NI}^{\ast}=1.1224{\hspace{0.17em}}_{-0.00003}^{+0.00007}$. This would correspond, for a room temperature nematic, to an uncertainty of ≈0.03 K, a significant improvement on previous estimates in literature. We have attempted a determination of the various critical exponents, which, even if roughly in agreement with experiments, turned out to be not sufficiently robust, as we discuss in §2. We estimated the divergence temperature of pretransitional ordering, finding ${T}_{C}^{\ast}\sim 1.1219$ from the heat capacity and from the correlation length obtained by an Ornstein–Zernike decay in the isotropic phases. We also evaluated the divergence of the pretransitional field susceptibility from the so-called second-order Kirkwood factor *g*_{2}, similarly to what is done in the analysis of real experimental data (e.g. [7,8,47]) by plotting ${T}^{\ast}/{g}_{2}$ versus ${T}^{\ast}$ and determining the intercept with the abscissa of its linear approximation. We found that the temperature dependence of ${T}^{\ast}/{g}_{2}$ (figure 5) shows a linear behaviour with different slopes for two temperature regions. The first, not too close to the NI transition $(1.1260\le {T}^{\ast}\le 1.1500$), is compatible with the predictions of the Landau phenomenological theory and gives an estimated ${T}_{c}^{\ast}\sim 1.1202$, while the second, close to the transition, where fluctuations become important, gives ${T}_{c}^{\ast}\sim 1.1219$, closer to ${T}_{NI}^{\ast}$. We have also discussed the shape and size of the ordered domains forming as *T*_{NI} is approached from above, determining the domain size distribution *p*(*s*), with an application of the Westin metric for the analysis of second rank tensors [44] to the ordering matrix [45,46]. Close to the NI transition, we find the distribution to be compatible with a power law, $p(s)\propto {s}^{\tau}$, with *τ* ∼ −1.4, over a broad range of domain sizes *s*. In summary, we believe that this new attack and fresh MC simulation data on the simplest model of the nematic-isotropic transition can help in clarifying some elements of this yet puzzling phase change [9] showing at the same time first- and second-order transition features.

### Data accessibility

This article has no additional data.

### Authors' contributions

G.S. and C.Z. designed the study and drafted the manuscript. G.S. carried out the simulations. Both authors read and approved the manuscript.

### Competing interests

We declare we have no competing interests.

### Funding

Slovenian Research Agency (grant nos P1-0099 and J1-9147).

## Acknowledgements

C.Z. would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality (under EPSRC grant no. EP/R014604/1) during the programme ‘Mathematical Design of New Materials’ when work on this paper was started.

### References

- 1.
Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E . 1953 Equation of state calculations by fast computing machines.**J. Chem. Phys.**, 1087–1092. (doi:10.1063/1.1699114) Crossref, Web of Science, Google Scholar**21** - 2.
Wang F, Landau DP . 2001 Efficient, multiple–range random walk algorithm to calculate the density of states.**Phys. Rev. Lett.**, 2050–2053. (doi:10.1103/PhysRevLett.86.2050) Crossref, PubMed, Web of Science, Google Scholar**86** - 3.
Jayasri D, Sastry VS, Murthy KP . 2005 Wang-Landau Monte Carlo simulation of isotropic-nematic transition in liquid crystals.**Phys. Rev. E**, 036702. (doi:10.1103/PhysRevE.72.036702) Crossref, Web of Science, Google Scholar**72** - 4.
de Gennes PG, Prost J . 1993**The physics of liquid crystals**. Oxford, UK: Oxford University Press. Google Scholar - 5.
Stinson TW, Litster JD . 1970 Pretransitional phenomena in the isotropic phase of a nematic liquid crystal.**Phys. Rev. Lett.**, 503–506. (doi:10.1103/PhysRevLett.25.503) Crossref, Web of Science, Google Scholar**25** - 6.
Ruessink BH, Barnhoorn J, Bulthuis J, Maclean C . 1988 Electric-field NMR of pretransitional effects in liquid crystals - a solute study.**Liq. Cryst.**, 31–41. (doi:10.1080/02678298808086347) Crossref, Web of Science, Google Scholar**3** - 7.
Attard GS, Beckmann PA, Emsley JW, Luckhurst GR, Turner DL . 1982 Pretransitional behavior in nematic liquid crystals. A nuclear magnetic resonance study.**Mol. Phys.**, 1125–1129. (doi:10.1080/00268978200100861) Crossref, Web of Science, Google Scholar**45** - 8.
Luckhurst GR . 1988 Pretransitional behaviour in liquid crystals.**J. Chem. Soc. Faraday Trans.**, 961–986. (doi:10.1039/f29888400961) Crossref, Google Scholar**84** - 9.
Mukherjee PK . 1998 The puzzle of the nematic-isotropic phase transition.**J. Phys. Cond. Matter**, 9191. (doi:10.1088/0953-8984/10/41/003) Crossref, Web of Science, Google Scholar**10** - 10.
Coles HJ, Jennings BR . 1978 Optical and electrical Kerr effect in 4-normal-pentyl-4’-cyanobiphenyl.**Mol. Phys.**, 1661–1673. (doi:10.1080/00268977800102661) Crossref, Web of Science, Google Scholar**36** - 11.
Keyes PH . 1978 Tricritical behavior at the isotropic-nematic transition.**Phys. Lett. A**, 132–134. (doi:10.1016/0375-9601(78)90043-9) Crossref, Web of Science, Google Scholar**67** - 12.
Blachnik N, Kneppe H, Schneider F . 2000 Cotton-Mouton constants and pretransitional phenomena in the isotropic phase of liquid crystals.**Liq. Cryst.**, 1219–1227. (doi:10.1080/02678290050122079) Crossref, Web of Science, Google Scholar**27** - 13.
Lebwohl PA, Lasher G . 1972 Nematic liquid crystal order. A Monte Carlo calculation.**Phys. Rev. A**, 426–429. (doi:10.1103/PhysRevA.6.426) Crossref, Web of Science, Google Scholar**6** - 14.
- 15.
Chen K, Ferrenberg AM, Landau DP . 1993 Static critical behavior of three-dimensional classical Heisenberg models: a high-resolution Monte Carlo study.**Phys. Rev. B**, 3249–3256. (doi:10.1103/PhysRevB.48.3249) Crossref, Web of Science, Google Scholar**48** - 16.
Holm C, Janke W . 1993 Critical exponents of the classical three-dimensional Heisenberg model: a single-cluster Monte Carlo study.**Phys. Rev. B**, 936. (doi:10.1103/PhysRevB.48.936) Crossref, Web of Science, Google Scholar**48** - 17.
Holm C, Janke W . 1997 Critical exponents of the classical Heisenberg ferromagnet.**Phys. Rev. Lett.**, 2265–2265. (doi:10.1103/PhysRevLett.78.2265) Crossref, Web of Science, Google Scholar**78** - 18.
Campostrini M, Hasenbusch M, Pelissetto A, Rossi P, Vicari E . 2002 Critical exponents and equation of state of the three-dimensional Heisenberg universality class.**Phys. Rev.****B**, 65. (doi:10.1103/PhysRevB.65.144520) Google Scholar - 19.
Luckhurst GR, Zannoni C . 1977 Why is the Maier-Saupe theory of nematic liquid crystals so successful?**Nature**, 412–414. (doi:10.1038/267412b0) Crossref, Web of Science, Google Scholar**267** - 20.
Challa MSS, Landau DP, Binder K . 1986 Finite-size effects at temperature-driven first-order transitions.**Phys. Rev. B**, 1841. (doi:10.1103/PhysRevB.34.1841) Crossref, Web of Science, Google Scholar**34** - 21.
Binder K . 1981 Static and dynamic critical phenomena of the two-dimensional q-state Potts model.**J. Stat. Phys.**, 69–86. (doi:10.1007/BF01007636) Crossref, Web of Science, Google Scholar**24** - 22.
Peczak P, Landau DP . 1988 Finite size effects at first-order phase transitions revisited. In*Computer Simulation Studies in Condensed Matter Physics: Recent Developments*(eds DP Landau, KK Mon, H-B Schüttler). Berlin, Germany: Springer. Google Scholar - 23.
Peczak P, Landau DP . 1989 Monte Carlo study of finite-size effects at a weakly 1st-order phase-transition.**Phys. Rev. B**, 11932–11942. (doi:10.1103/PhysRevB.39.11932) Crossref, Web of Science, Google Scholar**39** - 24.
Binder K . 1979**Monte Carlo Methods in Statistical Physics**. Berlin, Germany: Springer. Crossref, Google Scholar - 25.
Binder K . 1987 Theory of first order phase transitions.**Rep. Progr. Phys.**, 783. (doi:10.1088/0034-4885/50/7/001) Crossref, Web of Science, Google Scholar**50** - 26.
Fabbri U, Zannoni C . 1986 Monte Carlo investigation of the Lebwohl–Lasher lattice model in the vicinity of its orientational phase transition.**Mol. Phys.**, 763–788. (doi:10.1080/00268978600101561) Crossref, Web of Science, Google Scholar**58** - 27.
Zannoni C . 1986 A Cluster Monte Carlo method for the simulation of anisotropic systems.**J. Chem. Phys.**, 424–433. (doi:10.1063/1.450155) Crossref, Web of Science, Google Scholar**84** - 28.
Zhang Z, Mouritsen OG, Zuckermann MJ . 1992 Weak first-order orientational transition in the Lebwohl-Lasher model for liquid crystals.**Phys. Rev. Lett.**, 2803–2806. (doi:10.1103/PhysRevLett.69.2803) Crossref, PubMed, Web of Science, Google Scholar**69** - 29.
Priezjev NV, Pelcovits RA . 2001 Cluster Monte Carlo simulations of the nematic-isotropic transition.**Phys. Rev. E**, 062702. (doi:10.1103/PhysRevE.63.062702) Crossref, Web of Science, Google Scholar**63** - 30.
Shekhar R, Whitmer JK, Malshe R, Moreno-Razo JA, Roberts TF, de Pablo JJ . 2012 Isotropic-nematic phase transition in the Lebwohl-Lasher model from density of states simulations.**J. Chem. Phys.**, 234503. (doi:10.1063/1.4722209) Crossref, PubMed, Web of Science, Google Scholar**136** - 31.
Allen MP, Tildesley DJ . 2017**Computer simulation of liquids**, 2nd edn. Oxford, UK: Oxford University Press. Crossref, Google Scholar - 32.
Vollmayr K, Reger JD, Scheucher M, Binder K . 1993 Finite size effects at thermally-driven first order phase transitions: a phenomenological theory of the order parameter distribution.**Zeitschrift für Physik B**, 113–125. (doi:10.1007/BF01316713) Crossref, Web of Science, Google Scholar**91** - 33.
Priezjev NV, Pelcovits RA . 2001 Disclination loop behavior near the nematic-isotropic transition.**Phys. Rev. E**, 031710. (doi:10.1103/PhysRevE.64.031710) Crossref, Web of Science, Google Scholar**64** - 34.
Cordoyiannis G, Pinto LFV, Godinho MH, Glorieux C, Thoen J . 2009 High-resolution calorimetric study of the phase transitions of tridecylcyanobiphenyl and tetradecylcyanobiphenyl liquid crystals.**Phase Trans.**, 280–289. (doi:10.1080/01411590902858720) Crossref, Web of Science, Google Scholar**82** - 35.
Barker JA, Watts RO . 1969 Structure of water; a Monte Carlo calculation.**Chem. Phys. Lett.**, 144–145. (doi:10.1016/0009-2614(69)80119-3) Crossref, Web of Science, Google Scholar**3** - 36.
Wolff U . 1989 Collective Monte Carlo updating for spin systems.**Phys. Rev. Lett.**, 361. (doi:10.1103/PhysRevLett.62.361) Crossref, PubMed, Web of Science, Google Scholar**62** - 37.
Eppenga R, Frenkel D . 1984 Monte-Carlo study of the isotropic and nematic phases of infinitely thin hard platelets.**Mol. Phys.**, 1303–1334. (doi:10.1080/00268978400101951) Crossref, Web of Science, Google Scholar**52** - 38.
Lee J, Kosterlitz JM . 1991 Finite-size scaling and Monte Carlo simulations of first-order phase transitions.**Phys. Rev. B**, 3265. (doi:10.1103/PhysRevB.43.3265) Crossref, Web of Science, Google Scholar**43** - 39.
Blackman D, Vigna S . 2019 Scrambled linear pseudorandom number generators. (http://arxiv.org/abs/1805.01407v2). Google Scholar - 40.
Press WH, Teukolsky SA, Vetterling WT, Flannery BP . 1997**Numerical Recipes in Fortran 77: The Art of Scientific Computing**, 2nd edn. Cambridge, UK: Cambridge University Press. Google Scholar - 41.
Thoen J . 1995 Thermal investigations of phase transitions in thermotropic liquid crystal. In*Liquid Crystals in the Nineties and Beyond*(ed. S Kumar), pp. 19–80. World Scientific, Singapore. Google Scholar - 42.
Stanley HE . 1971**Introduction to phase transitions and critical phenomena**. Oxford, UK: Oxford University Press. Google Scholar - 43.
Humphries RL, Luckhurst GR . 1982 Computer simulation studies of anisotropic systems. VII. Pretransitional effects in nematogens.**Proc. R. Soc. Lond. A**, 307–318. (doi:10.1098/rspa.1982.0103) Link, Web of Science, Google Scholar**382** - 44.
Westin C-F, Maier SE, Mamata H, Nabavi A, Jolesz FA, Kikinis R . 2002 Processing and visualization for diffusion tensor MRI.**Med. Image Anal.**, 93–108. (doi:10.1016/S1361-8415(02)00053-1) Crossref, PubMed, Web of Science, Google Scholar**6** - 45.
Callan-Jones AC, Pelcovits RA, Slavin VA, Zhang S, Laidlaw DH, Loriot GB . 2006 Simulation and visualization of topological defects in nematic liquid crystals.**Phys. Rev. E**, 061701. (doi:10.1103/PhysRevE.74.061701) Crossref, Web of Science, Google Scholar**74** - 46.
Skačej G, Zannoni C . 2014 Molecular simulations shed light on supersoft elasticity in polydomain liquid crystal elastomers.**Macromolecules**, 8824–8832. (doi:10.1021/ma501836j) Crossref, Web of Science, Google Scholar**47** - 47.
Attard GS, Emsley JW, Khoo SK, Luckhurst GR . 1984 Pre-transitional behavior in the isotropic phase of nematogenic mixtures: an NMR investigation.**Chem. Phys. Lett.**, 244–248. (doi:10.1016/0009-2614(84)85022-8) Crossref, Web of Science, Google Scholar**105**