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

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₂〉. 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′-n-octyl-4-cyano biphenyl (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 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_BT/ϵ = 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

where ϵ_ij is a characteristic interaction energy equal to a positive constant, ϵ, when i, j are nearest neigbours and 0 otherwise, while P₂(x) ≡ (3x² − 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³ 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³ 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⁻⁴ in reduced temperature scale, $T^{\ast }=k_{B}T/ϵ$, 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_BT). 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ϵ)}^{-1}\sum _{⟨i,j⟩}{U}_{ij}$ and its average over 1.2 × 10⁶ MC production cycles (1 cycle corresponding to N attempted individual MC moves) is plotted versus reduced temperature $T^{\ast }\equiv k_{B}T/ϵ$ in figure 1a, with a close-up of the transition region in figure 1b. 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 1b). 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.

• Download figure
• Open in new tab
• Download PowerPoint

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₂〉 shown in figure 2a,b. In the isotropic phase one has 〈P₂〉 ∼ 0, while a non-zero 〈P₂〉 is a signature of nematic order. Note that the above algorithm always gives a positive contribution to the 〈P₂〉 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 $⟨P_2^m⟩$, also shown in figure 2a,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.

• Download figure
• Open in new tab
• Download PowerPoint

The order parameters 〈P₂〉, $⟨P_2^m⟩$ show a discontinuous jump at the NI transition. A more detailed look at the transition region (figure 2b) reveals, like for the energy in figure 1b, 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[⟨{(U^{\ast })}^2⟩-{⟨U^{\ast }⟩}^2]/{(T^{\ast })}^2$. Alternatively, it can also be calculated directly from its definition, $c_V^{\ast }=\text{d}⟨U^{\ast }⟩/\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 3a,b, exhibit a pronounced $c_V^{\ast }$ peak in the transition region.

• Download figure
• Open in new tab
• Download PowerPoint

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⁶ MC production cycles is observed. While for the relatively small system sizes of N = 30³ 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³ considered in [29] and is to some extent effective even for some N = 100³ 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³. 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 $⟨U^{\ast }⟩$ from simulation, and integrate the relation $⟨U^{\ast }⟩=\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,

starting the calculation from the available low- and high-temperature analytic results. Then, in principle, an intersection of the two curves would give the desired 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 $⟨U^{\ast }⟩({\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³ sample, within four-digit precision, as $T_{NI}^{\ast }={1.1224}_{-0.00003}^{+0.00007}$. Moreover, from figures 1b–2b, 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/ϵ\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₂〉 in the nematic phase, see figure 2a,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 3a,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₂〉), 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₂〉, 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₀(T) is a reference free energy that depends only on temperature. In this framework, the NI transition takes place at $T_{NI}=T_{\ast }+2b^2/9ac$, while the superheating temperature limit of the nematic phase is reached at $T^{†}=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^{†}-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₂(r) = 〈P₂(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₂(r). At each temperature, the resulting G₂(r) were fitted to the Ornstein–Zernike form G₂(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.

• Download figure
• Open in new tab
• Download PowerPoint

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₂ [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₂ can be calculated directly from the simulated G₂(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₂(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 5a,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 5b).

• Download figure
• Open in new tab
• Download PowerPoint

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

• Download figure
• Open in new tab
• Download PowerPoint

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₂(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⁴, 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 }$.

• Download figure
• Open in new tab
• Download PowerPoint

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³ and predicted $T_{NI}^{\ast }=1.1232\pm 0.0001$, which agrees with the early analysis by Fabbri and Zannoni in an N = 30³ sample [26]. Furthermore, Shekhar et al. [30] considered lattices up to N = 40³ 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³, 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³) 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.

• Download figure
• Open in new tab
• Download PowerPoint

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