Phase behaviour of the confined lattice gas Lebwohl-Lasher model

The phase behaviour of the Lebwohl-Lasher lattice gas model (one of the simplest representations of a nematogenic fluid) confined in a slab is investigated by means of extensive Monte Carlo simulations. The model is known to yield a first order gas-liquid transition in both the 2D and 3D limits, that is coupled with an orientational order-disorder transition. This latter transition happens to be first order in the 3D limit and it shares some characteristic features with the continuous defect mediated Berezinskii-Kosterlitz-Thouless transition in 2D. In this work we will analyze in detail the behaviour of this system taking full advantage of the lattice nature of the model and the particular symmetry of the interaction potential, which allows for the use of efficient cluster algorithms.


Introduction
The Lebwohl-Lasher (LL) model [1] is known to be one of the simplest systems that can reproduce the isotropic-nematic transition, which is a key feature in the physics of liquid crystals, ubiquitous materials in today's technology. It is actually nothing but a lattice version of the somewhat older continuum Maier-Saupe model [2]. Prof. Myroslav Holovko, whose 70th birthday we are celebrating with this Festschrift, was one of the first to study the solution of anisotropic integral equation approaches to study the Maier-Saupe fluid and its order-disorder transitions [3], and he has recently published a study in which the system is considered in terms of a field theoretical approach [4]. In this work, as a tribute to his numerous contributions to the field of statistical mechanics of fluids and phase transitions, we will present a computer simulation study of the phase behaviour of a Lebwohl-Lasher lattice (LLL) gas under confinement. This model is an extension of the previous work by us [5] in which the nature of the orientational transitions of the LL model was analyzed in depth, and its close connection with the defect mediated Berenzinskii-Kosterlitz-Thouless (BKT) transition [6,7] was investigated. Here, we add translational degrees of freedom in terms of a lattice occupation variable. It is worth to recall that the precise nature of the apparently ordered phase in the two dimensional version of the LL model has been very much discussed (see for instance reference [8] and references therein). In this paper we adopt the view of Reference [5], in the sense that these systems exhibit a transition between an isotropic and a quasi-nematic phase endowed with quasi-long range orientational order. We will confirm the findings for the confined continuum Maier-Saupe fluid which we investigated using computer simulations [9] and extend the calculations to much larger system sizes that are computationally unfeasible in continuum models. The coupling of the first order gas-liquid transition induced by the net attraction between the spins and the BKT-like transition of the confined fluid will be illustrated for various degrees of confinement. This BKTlike orientational transition evolves gradually into a weakly first order transition in the 3D bulk limit [10].
The rest of the paper is organized as follows. In the next section we describe the model and the Monte Carlo methods used to study the system. The last section is devoted to a presentation of the most significant results and conclusions.

Model and methods
Our LL model is the simplest version of a nematogenic fluid, in which every molecule interacts solely with molecules placed at neighbouring sites of the lattice. Every lattice site can be either occupied or empty, and this is controlled by an occupation number variable, n i , which can take the value 1 or 0 depending on whether there is a molecule at the site i or not. The total potential energy of the system can thus be expressed as U = −ε 〈i j 〉 n i n j P 2 (s i s j ), (2.1) where ε is the coupling parameter that defines the energy scale (ε > 0) and the reduced temperature T * = k B T /ε (where k B is Boltzmann's constant as usual), s i and s j are unit vectors that describe the orientation of molecules i and j , P 2 is the second degree Legendre polynomial and finally 〈i j 〉 indicates that the summation is restricted to nearest neighbour (NN) pairs of sites. Our model consists in a slab of cubic lattice, with periodic boundary conditions in the x, y directions and a width of H sites along the z direction as in the case of reference [5], where all sites were occupied. The model thus will have a total of L × L × H sites. According to the previous results [11], we expect to find a first-order (liquid-vapor-like) transition at low temperatures, and a continuous transition between an isotropic and a quasi-nematic phase at higher temperature. In the range of temperatures where the transition is continuous we have performed, for each size H , simulations for a series of L values in order perform a finite size scaling analysis, namely L = 10, 20, 30, 40, 60, 80, 90, 100.
As in the previous works, for continuous transitions we have performed Monte Carlo simulations using a combination of a local update algorithm [11,12] and cluster moves [11,13,14] so as to minimize critical slowdown effects when approaching the critical temperature. The temperature range in which first-order transitions are expected to occur was analysed by means of simulations using Grand-Canonical Wang-Landau (GCWL) methodology, following the prescription of reference [9] for continuum models. We refer the reader to our previous works [5,9,11,15] for technical details about the precise implementation of the methods. Once more, the isotropic-quasi-nematic transition is monitored by means of the largest eigenvalue, λ + , of Saupe's tensor [16] where N is the number of spins, and α and β refer to the x, y, z components of the unit vector s i . For a fixed value of the chemical potential, µ, we can define a size dependent pseudocritical temperature T c (µ, L, H ), in terms of the fluctuation of the order parameter Using the maxima of the L-dependent susceptibility as the estimate of T c (µ, L), as in previous works, we have assumed a BKT-like finite size scaling of the form [17,18]    In the case of the first order gas-liquid transition, we have carried out computer simulations for different temperatures and systems sizes using the GCWL method as follows: • For a given T and L, the Helmholtz energy is calculated for a series of densities in the fluid regime.
• The possible phase transition is located by determining the value of the chemical potential µ e (L, T ) that maximizes the density fluctuations. Paying attention to the grand canonical density distribution functions [11] P (ρ; µ e , T, L), and in particular to the dependence of their shape with L, it is possible to probe the existence of the first order transition at the corresponding temperature, and eventually to estimate the apparent coexistence densities of the vapour and liquid phases, ρ v (L, T ), and ρ l (L, T ). These system-size dependent densities can be used to extrapolate the results to the thermodynamic limit (L → ∞).
Finally, the results for µ tc (L), T tc (L) and ρ tc (L) are fitted to second degree polynomials of 1/L, where the linear term is taken from the analysis of the related planar Maier-Saupe model [11] and the quadratic term accounts for the effects of confinement [9]. In this way we can finally obtain the estimates for L → ∞.

Results
As mentioned before, calculations have been carried out for a series of system sizes and values of H = 1, 2, 3, 4, 8 and bulk 3D system. The system can be thought of as either a slab in vacuum or a fluid confined in an inert slit pore (wall-particle interactions reduced to hard core exclusion). In figure 1 we plot the order parameter and susceptibility results for the limiting case H = 1 (a 2D LL fluid). As it is characteristic of BKT-like transitions [5], the order parameter, although exhibits a clear jump across a transition temperature, tends to vanish with an increasing sample size. On the other hand, the susceptibility shows an apparent divergence, both at the transition temperature and at temperatures below. The values of the fraction of percolating clusters, X , shown in the upper graph of figure 2 exhibit a clear jump at the transition temperature. The corresponding grand canonical constant volume heat capacity plotted in the lower graph of the same figure (where u = U /V is the potential energy density), exhibits no size dependent divergence, but most probably tends to build up a cusp singularity characteristic of BKT-like transitions [17]. Again, similar results are obtained for temperatures/chemical potentials above the transition temperature/chemical potential when analyzing the results for other H values, until in the limit of H → ∞ one recovers the bulk behaviour, i.e., a first order phase transition as shown in reference [5]. For finite pore widths below a certain multicritical temperature T tc , one encounters first order transitions. This is illustrated in figure 3 and in figure 4, where one can appreciate large variations of the density and the energy with the chemical potential correlated with the BKT-like transition for a temperature slightly above T tc . These multicritical temperatures (the highest T at which the isotropic-quasinematic transition is coupled with a first order gas-liquid transition) are estimated using the approach indicated in the previous section, leading to the results of table 1.
The T −ρ phase diagram for the LLL model under varying degrees of confinement is presented in figure 5. We see how the relatively flat (Ising 2D) behaviour of the vapour-liquid curve for H = 1 gradually evolves into the shape of the bulk 3D LL isotropic-nematic transition in which the orientational transition is fully coupled to the vapour-liquid transition for all temperatures. If one tries to correlate the value of the multicritical temperatures with the pore, we again find that the behaviour deviates somewhat from the Kelvinlike scaling T c (bulk) − T c (H ) ∝ 1/H as it was the case in the confined Maier-Saupe continuum case. Due to the presence of the order-disorder transition, the multicritical points cannot be directly mapped into vapour-liquid critical points, in particular, as the width of the pore increases. As a matter of fact, the H = 1 phase diagram might also appear to be compatible with the presence of a critical endpoint, instead of a tricritical point, but on the basis of the behaviour for H > 1 one can think that most likely the BKT line for H = 1 (solid red circles) bends backwards as T is lowered to hit the vapour-liquid equilibrium curve right at the top (the VL critical point). As H is increased this behaviour is more apparent.
In summary, we have presented a detailed analysis of the phase behaviour of the confined Lebwohl-Lasher lattice model. We have shown how the confinement transforms the first order isotropic-nematic transition of the bulk 3D LL system into a continuum BKT transition between an isotropic phase and a phase with local orientational order (quasinematic), which develops into a first order vapour-liquid transition below a given multicritical temperature. These findings confirm and extend those for the continuum Maier-Saupe model with a more reliable finite size scaling analysis.