Nematic-isotropic transition in a weakly diluted lattice model: Monte Carlo study

We consider a weakly diluted lattice model of elongated particles undergoing a rst order nematic-isotropic transition in a pure bulk case. The lattice of 20 3 sites with a 5% dilution is studied by Monte Carlo simulation and by applying the Ferrenberg-Swendsen histogram technique. The shift of the transition temperature and decrease of the latent heat and the order parameter due to dilution are calculated. The heat capacity and the susceptibility peaks are essentially suppressed as compared to a pure case. The comparison to the experiments on liquid crystals connned in silica aerogels and porous glasses is done.


Introduction
The study of the nematic-isotropic (NI) transition in liquid crystals (LCs) in the presence of impurities has several reasons.In general, the experimental systems studied are prepared with some level of purity and thus, have impurities of different kind which lead to rounding of the NI transition.As a result, slightly broadened peaks for the heat capacity as well as rounded edges in the enthalpy jumps near the transition are observed [1].In this context the influence of dilution on the first order transitions in three dimensions was examined theoretically by Imry and Wortis [2].They showed that quenched impurities can round the transition, where the strength of the rounding depends on the fluctuation in the impurity density and the interfacial energy between the phases.Berker [3] showed that discontinuities can be suppressed in all dimensions, dependent on the symmetry of the system and the strength of the random field.These results have been reexamined recently by Cardy [4] on a base of the q-state Potts model with quenched random impurities.Previous studies of the first order transitions with the different type of quenched disorder by computer simulation compared the shift of the transition temperature and the rounding of the c Ja.M.Ilnytskyi specific heat with the pure case [5][6][7][8][9].It is found that the transition can change from first to second order above some nonzero threshold value of impurity concentration [6,8].
The other reason to study the effect of dilution is to probe it as a simple description of the porous confining media.The main phenomena observed in the experiments on LCs in a porous media, namely the lowering of the transition temperature, the decrease of the transition enthalpy, suppression and broadening of the heat capacity peak, suppression of the orientational order [10][11][12][13][14][15] recall the effect of dilution.However, the strength of these effects depends crucially on the porosity of the confinement material, its structure and the surface treatment of its pores.Being a simple model, the dilution neglects some of these aspects, but retains two important ones, such as finite-volume effect and a connectivity of the pores.
The aim of this report is to perform a Monte Carlo (MC) study of the NI transition in a lattice model proposed earlier [16] in the presence of a quenched random dilution.The particles in this model interact via the angular part of the Berne-Pechukas (BP) potential [17] derived from the overlap integral of two ellipsoidal Gaussians of a certain elongation a.The elongation of the particles considered in this report is a = 3.In this case reasonable values for the latent heat and the order parameter at the NI transition in pure bulk systems are reported [16].We restrict ourselves to a single lattice size of 20 3 sites with a 5% dilution.This situation may be addressed to a LC in a high porosity random confining media.The present study has two main objectives, namely, to find out how dilution affects the thermodynamics near the NI transition, whether the transition remains first order, and how these changes compare with experiments on LCs in high porosity confining media.To improve the accuracy of our analysis the Ferrenberg-Swendsen (FS) histogram technique [18] is used.

Monte Carlo simulations
We consider a lattice model with N = 20 3 sites.A part of these sites N m = cN are free and are concerned as impurities.Each of the other N f = N(1 − c) sites represents a particle with certain orientation given by a unit vector ûi .The nearest neighbours interaction between the particles has the form [16] V where a is the elongation parameter (we use a = 3 in the present study), ǫ is the energy scale, χ = (a 2 − 1)/(a 2 + 1) is the anisotropy parameter, and cos θ ij = ûi ûj .This potential is derived from the overlap of two ellipsoidal Gaussians with the elongation a [17] and can be considered as some kind of mapping between hardcore and continuous models [19].It is interesting to note that the potential (1) for a = 3 perfectly coincides with the potential derived from the excluded volume of two spherocylinders [19] and expanded up to P 6 (cos θ ij ) term [20].In the limit of small anisotropy χ ≪ 1 the exact form of the Lebwohl-Lasher potential [21] is reproduced.As a increases the potential (1) becomes more anisotropic (see [16]).It has to be noted that a similar effect can be achieved by considering higher order P 2n (cos θ ij ) terms [22][23][24] in addition to P 2 (cos θ ij ) in the Lebwohl-Lasher potential.
A standard Metropolis algorithm is used to simulate this model in the vicinity of the NI transition.The orientation of each particle ûi is changed by adding a vector l of random orientation and controlled lenght [25] and then normalizing û ′ i = ûi + l back to unity.The lenght of l is adjusted automatically during simulations to achieve an acceptance ratio of about 0.4.The dimensionless temperature T * = k B T /ǫ is used in our simulation, and the dimensionless single-particle internal energy U * and the order parameter (where θ i is the angle between the long axis of the i-th particle and a director) are evaluated after each simulational cycle.The method of Viellard-Baron [26] is used to calculate S. Firstly, the short scanning runs (up to 10 5 MC cycles) were performed along the entire interval of temperatures including the NI transition point.The more extended run of 5 • 10 5 MC cycles was done then for the temperature at which the best coexistence of two phases had been observed.The histograms of the energy and the order parameter distribution are built and the FS reweighting [18] is used to locate the NI transition point in a first approximation.This temperature is then used for the final extended run of minimum 10 6 MC cycles which gives the final histograms and the quantities of interest.We calculate the dimensionless heat capacity C * v = C v /(k B N f ) and the susceptibility χ * = χǫ/N f defined via the fluctuational formulae: Binder's fourth cumulant [27] of the energy fluctuations is also calculated since it can be used in the investigation of the order of the transition.The method is the same both in a pure and a diluted cases.Three different configurations of 5% quenched dilutions are used in these simulations and the averaging over them is performed.
Due to finite-size effects [28] the peaks in the heat capacity and the susceptibility are broadened and have maxima at slightly different temperatures T * 1,NI and T *

2,NI
respectively.Binder's fourth cumulant has a minimum at T * 3,NI .The histograms of the energy distribution have double maxima of equal height at a certain temperature T * 4,NI being close to all temperatures mentioned above.In the limit of an infinite system, temperatures of all these estimates coincide (see, for example, [29]) giving the transition temperature T * NI with high accuracy.In our case of a single lattice size we are more interested in the shift of the transition temperature due to dilution than in its exact value.We obtain substantial shifts for all four estimates T * 1−4,NI caused by a dilution.However, the ratio turned out to be the same for all m = 1, 2, 3, 4 within the accuracy of our calculations.
The maxima for C * v max and χ * max are suppressed by dilution (see figure 1).We obtain the ratios To obtain the discontinuous latent heat of the NI transition ∆H NI consider the energy distribution histogram at temperature T * 4,NI .It is known that the energy distribution of a system close to a first order transition is approximated reasonably well by a double Gaussian [30], however, we obtain much better fits using double non-Gaussians of the form: where the expected values in the nematic and the isotropic phase U nem and U iso and the fitting coefficients are obtained numerically using the least-squares method (the accuracy χ 2 is of order 0.01).So, the dimensionless one-particle latent heat was then estimated as ∆H * NI = U * iso − U * nem (see figure 2).
The same method is used to estimate the order parameter at the transition S NI .In this case the order parameter distribution is considered at T * 2,NI .At this temperature the maximum correspondent to the isotropic phase is very low and it is fitted by a Gaussian, and the one correspondent to the nematic phase by a non-Gaussian: (see figure 2) and we assumed that S NI = S nem .
-0.8 -0.7 -0.6 -0.5 -0.4 -0.We obtain the following ratios for the latent heat and the order parameter at the transition in a pure and a diluted model: Thus, the first order NI transition which takes place in the pure version of the model considered here is essentially smeared due to the presence of a weak 5% dilution.The transition temperature is shifted, the maxima for the heat capacity and the susceptibility are suppressed.The latent heat and the order parameter are essentially lower as in the corresponding pure case.However, there are indications that the transition retains its first order nature.To prove this unambiguously the finite size scaling analysis has to be performed [31].

Discussion
The influence of a dilution on the thermodynamics near the NI transition can be compared with the corresponding experiments on LCs confined in a high porosity media.Particularly, Wu et al. [12] have studied the NI transition in 8CB LC confined to silica aerogels with a different porosity.For ρ = 0.08 g cm −3 aerogel density (which corresponds roughly to the 5% volume fraction of impurities in our model) a shift of T NI of -0.45 • C was observed, thus giving T gel NI /T pure NI = 0.9986.The shift obtained in our simulations and given by a ratio (3) is much more pronounced.We can also compare the suppression of the heat capacity maxima by progressive increasing of the aerogel density ρ.As follows from the experiments [12], the excess heat capacity ∆C p max decreases approximately linearly with an increase of ρ for ρ < 0.36 g cm −3 .For the aerogel densities ρ 1 = 0.08 g cm −3 and ρ 2 = 0.17 g cm −3 we obtain the ratio ∆C p max (ρ 2 )/∆C p max (ρ 1 ) = 0.51, which is close to the ratio 0.45 (4) obtained in our simulations.
The other experiments on the NI transition in 8CB LC confined to porous glasses were done by Iannacchione et al. [15].In the case of macroporous confinement (1000 Å mean pore size) a shift of the transition temperature of -2.05 • C was observed.This gives a ratio T glass NI /T pure NI = 0.993 which is again higher than the ratio found in (3).The latent heat ∆H glass is reduced as compared to a pure case and ∆H glass /∆H pur = 0.74, the suppression of the heat capacity gives a ratio ∆C p max (glass)/∆C p max (pure) = 0.65.These values are higher than the ones obtained in our simulations (7,4).This, probably, is due to the fact that a macroporous confinement with 1000 Å mean pore size corresponds to a much weaker dilution than 5% considered in the simulation.
We obtain a suppression of the first order NI transition in a weakly diluted lattice model of elongated particles as compared to the pure case.The effects observed in the simulations agree qualitatively well with the experimentally observed ones in silica aerogels and support previous computer simulations of the other models [6,8].However, the shift of the transition temperature is essentially overestimated.This can be explained by the fact that a site in the lattice model describes a group of the real molecules rather than a single one.Then the dilution at one site will destroy 6 bonds for a sc lattice and the energy of 6 surrounding particles (groups of molecules) become essentially underestimated, so the transition temperature will be shifted too much.In this context we guess that a 5% dilution corresponds effectively to much higher density aerogel than ρ = 0.08 g cm −3

Figure 1 .
Figure 1.Suppression of the heat capacity C * v and the susceptibility χ * peaks in the vicinity of the NI transition due to a weak dilution.

Figure 2 .
Figure 2. Histograms of the energy distribution at T * 4,NI (on left) and the order parameter distribution at T * 2,NI (on right) and their fits by the forms (5,6).The expected values U * nem , U * iso and S nem are shown.