Nematic phase transitions in two-dimensional systems

Simulations of nematic-isotropic transition of liquid crystals in two dimensions are performed using an O(2) vector model characterised by non linear nearest neighbour spin interaction governed by the fourth Legendre polynomial $P\_4$. The system is studied through standard Finite-Size Scaling and conformal rescaling of density profiles or correlation functions. The low temperature limit is discussed in the spin wave approximation and confirms the numerical results, while the value of the correlation function exponent at the deconfining transition seems controversial.


Ordering in two dimensions
In the context of phase transitions, two-dimensional models exhibit a very rich variety of typical behaviours, ranging from conventional temperature-driven second order phase transitions (e.g. Ising model) to first-order ones (e.g. q > 4−state Potts model), with also specific properties of models with continuous global symmetry which may present defect-mediated topological phase transitions (e.g. XY model) or even no transition at all (e.g. Heisenberg model). Models of nematic-isotropic orientational phase transitions belong to this latter category of systems displaying a continuous symmetry. Ordering in low dimensional systems is likely to be frustrated by the strentgh of fluctuations. On qualitative grounds, consider e.g. how fluctuations develop within the framework of Landau theory when the temperature decreases from the high temperature paramagnetic phase toward the transition temperature. The response to a localized magnetic field applied at the origin, hδ(x), follows from Ginzburg-Landau functional minimization of the free energy F [m(x)] = dx( 1 2 am 2 (x) + 1 4 bm 4 (x) + K|∇m(x)| 2 − hδ(x)m(x)). After linearization and Fourier transform, it yields where ξ −2 = a/2K is the correlation length. The responsem(k) is here proportional to the correlation function Fourier transformG(k) and the fluctuations are measured through The latter sum is converted to an integral from 0 to some cutoff Λ. It is diverging with ξ in 1d ( kG (k) ∼ ξArctan (Λξ)) and in 2d ( kG (k) ∼ ln(Λξ)) while it is bounded in 3d ( kG (k) ∼ const.). The divergence in 1d prevents any long range ordering, while stable ordered state is not forbidden in three-dimensional systems.
In the intermediate case, due to the logarithmic diverging behaviour of the intergral it is less obvious to conclude and more refined analysis is required. In his famous book on phase transitions, Cardy uses a simplified version of Peierls argument on the existence of a phase transition at finite temperature in 2d in the case of discrete symmetry and extends the argument to continuous symmetry, showing that ordered ground state is unstable with respect to thermal fluctuations in this latter situation. Consider a spin system with only nearest-neighbour interactions −JS i S j and assume that the spins are represented by classical n−component vectors. An ordered ground state may be stabilized e.g. by symmetry breaking fields at some boundaries of the system. The variation of internal energy when a droplet of typical size l with spins progressively tilted in such a way that at the center of the droplet the spins are pointing opposite to the direction of the field is of order of O(l 2 ) × J|S| 2 (π/l) 2 where π/l is the nearest neighbour spin disorientation. This result follows from integration over the droplet volume O(l 2 ). Considering that the entropy is measured by the number of possible closed loops of size O(l) in 2d, the entropy of the droplet is estimated as k B ln µ l where µ < z − 1 (z is the coordination of the lattice), so that eventually the free energy variation is ∆F (l) = π 2 J|S| 2 − k B T l ln µ. At any non zero temperature, increasing the size l of the droplet stabilizes the system and a spontaneously ordered ground state is thus impossible. In the case of discrete symmetry, e.g. Ising model, the energy balance would be associated to the interface only and we would get ∆F (l) = 2J|S| 2 l − k B T l ln µ, showing that a transition to a stable ordered ground state is expected at a temperature around The two examples are illustrated in Fig. 1. In the latter situation, a conventional Ð Ð Figure 1. Evaluation of the free energy of a disordering droplet in two dimensions in the case of continuous symmetry (left) and discrete symmetry (right).
phase transition toward an ordered phase is expected at finite temperature, while such an ordered phase may only be encountered at zero temperature in the first example [1,2]. On the other hand, an unconventional phase transition toward a quasi-long-range ordered state may take place at finite temperature, as we discuss in the next section.

Two-dimensional electrodynamics and the XY model
The celebrated XY model, admits a phenomenological description in terms of Coulomb gas. This is a wellknown description, first given by Berezinskiȋ and Kosterlitz and Thouless [3][4][5] and it is worth reminding its essential steps here. Before considering this model, imagine a point-charge q located at site r 0 in a two-dimensional space, ρ(r) = qδ(r − r 0 ). The corresponding Coulomb potential, solution of Poisson equation ∇ 2 φ(r) = −ε −1 0 ρ(r) might be written where a is chosen such that φ(a) = 0. Consider now a spin system living in twodimensional space and let call u = ∇θ the distorsion field, where θ(r) is the phase field defined by S(r) = (cos θ(r), sin θ(r)). Due to the periodicity of the phase field, the distorsion field should obey the following relation C(r 0 ) u dl = 2π × integer where the contour integral is taken along a counterclockwise closed path around the point r 0 . Using Stokes theorem, we may also write The general solution of this equation consists of two terms, u = ∇ψ − ∇ × (ẑφ). The curl of the first term being identically zero, we are led to a Poisson equation for the singular term ∇ 2 φ = 2π × integer × δ(r − r 0 ) where the integer ≡ n plays the role of the total (topological) charge enclosed by the contour C. From an energetic point of view, we may consider the kinetic energy which, up to an inessential constant, is the continuous approximation of equation (3) with K = βJ. After decomposition in ψ and φ parts, the cross-term vanishes and we get two independent contributions, spin waves defects (8) This is a fundamental relation which is the basis of the factorization of the partition function in a spin wave contribution and a defect (vortex) contribution, Z = Tr exp(−βH) = Z SW Z V . In the lattice version, the spin wave contribution to the Hamiltonian βH SW = 1 2 K (r,r ′ ) (ψ r − ψ r ′ ) 2 is quadratic in Fourier space, (µ is the lattice spacing and ψ q the Fourier component of ψ r ). This expression leads to the Gaussian model which implies [6] that where a temperature-dependent spin-spin critical exponent is found, η SW = (2πK) −1 .
The low temperature (LT) phase of the XY model is a quasi-long-range ordered phase (QLRO), or a critical phase. When the temperature increases, vortices bounded in pairs produce further disordering of the system and a faster than linear increase of the exponent η with temperature up to the temperature where the transition takes place. The role of vortices is understood perturbatively through the calculation of the effective screened interaction energy between the topological charges. We consider a neutral Coulomb gas. Omitting a diverging contribution O(ln R/a) × i n i which would occur otherwise (i.e. for a non-neutral Coulomb gas), the vortices contribution to the Hamiltonian reads as where V (r i −r j ) is the Coulomb interaction between the charges. In the perturbative approach, we consider the effective interaction between charges n i and n j = −n i in the presence of another screening dipole. The following result comes out [7] where y 0 is the fugacity of a charge. The correction term y 2 0 dr r 3−2πK diverges at the deconfining transition 2πK c = 4 where the pairs of charges break. The presence of vortices increases the disordering of the system which, below T c , flows under renormalization to the zero-fugacity limit where the spin-wave limit is recovered with a renormalized temperature [5]. At the transition, the correlation function exponent takes a universal value The Heisenberg model (O (3)) on the other hand has no transition at any finite temperature (asymptotic freedom). An intuitive argument called escape to the third dimension is often reported. This is the observation than vortices cannot be stable for n ≥ 3, since an infinitesimal amount of energy is able to produce a spin-wavetype excitation which eliminates the localized defect. Also a renormalization scheme was proposed by Polyakov to treat the non linear σ−model close to the lower critical dimension d = 2 [8][9][10][11]. The flow of the coupling constant g (the temperature) under a change of scale s is governed by the β−function, β(g, d) = s ∂g which shows that the model is disordered at any temperature when n > 2 in two dimensions, since the coupling always decreases and flows to zero under renormalization.

Definition of the model and of the observables
We now come back to liquid crystals, the molecules of which may be idealized as long neutral rigid rods. They are likely to interact through electrostatic interactions, therefore Legendre polynomials appear for the description of the orientational transition between a disordered isotropic high temperature phase and an ordered nematic phase. Lattice models of nematic-isotropic transitions capture the essentials of this extremely simplified description. The molecules are represented by n−component unit vectors S r (also called "spins"), located here on the sites r of a square lattice. The interaction between molecules is restricted to the nearest neighbour pairs (r, r ′ ), the radial dependence being kept constant and the angular dependence entering through a k−th order Legendre polynomial 1 , P k (S r · S r ′ ), in terms of the scalar product between S r and S r ′ , S r ·S r ′ . A coupling parameter J measures the interaction intensity. One obtains the following Hamiltonian of a lattice liquid crystal, When the value of k in Eq. (14) is varied, new features may be expected, as in the case of symmetry-breaking magnetic fields h k cos kθ added to the XY model which change the phase diagram as investigated by José et al [12,7]. When the polynomial order k increases, one may expect a qualitative change in the nature of the transition, like in the case of discrete spin symmetries (Potts model) [13,14]. The value k = 2 was intensively studied. It still corresponds to the XY model for O(2) spin symmetry, while it leads to the RP 2 or Lebwohl-Lasher model [15] for 3−component spin vectors. The nature of the transition in this latter case is still under discussion: an early study of Kunz and Zumbach [16] reported numerical evidences in favour of a topological transition, but more recently, several authors argued in favour of the absence of any finite-temperature phase transition [17,18], like in the Heisenberg case. Our own previous contributions [19,20] support the first scenario with QLRO at low temperature like in the XY model, but one cannot exclude a finite -but extremely large -correlation length which exceeds the maximum size available in numerical simulations. A recent study reported extremely convincing new evidences in favour of a topological transition [21], the transition being driven by topologically stable point defects known as 1 2 disclination points. Eventually, in the large−n limit, there is a proof of asymptotic freedom for values of k (in the interaction term (1 + cos θ) k ) which do not exceed a critical k c ≃ 4.537... [18]. Above this value the transition becomes of first order, a result which does not violate Mermin-Wagner-Hohenberg theorem, since the correlation length is finite at the transition. For finite value of n, the question of the nature of the transition at high k is still a challenging problem, although there is a rigourous proof that the transition becomes of first-order for large enough values of k for arbitrary n ≥ 2 [22,23]. This observation suggests to inspect the effect of a higher value of k also for the O(2) model. In the context of orientational transitions in liquid crystals, Legendre polynomials rather than cos k θ interactions are introduced, and we are led to the Hamiltonian of Eq. (14).
In this report, we consider the behaviour of an Abelian spin model (O(2) rotation group) with P 4 −like spin interactions. We will refer to it as the P 4 O(2) model for simplicity. For 2−component vectors in a disordered phase, cos 2 θ = 1 2 and cos 4 θ = 3 8 . In order to keep the same symmetry in the interaction than in the P 4 O(3) model already considered in the literature [24,25], but to normalise it between 0 and 1 in the limits of completely disordered and completely ordered phases respectively, we modify slightly the Hamiltonian, considering pair interactions of the form Q 4 (x) ≡ AP 4 (x)+const = 8 55 (35x 4 −30x 2 + 15 8 ). The corresponding Hamiltonian is thus defined by with S r = (S x r , S y r ), |S r | = 1. A qualitative description of the transition is provided by the temperature behaviour of the energy density, the specific heat, the order parameter and the susceptibility. The internal energy is defined from the thermal average of the Hamiltonian density, u(T ) = (dL d ) −1 H and the specific heat follows from fluctuation dissipation theorem, C v (T ) = (L d T 2 ) −1 ( H 2 − H 2 ). Brackets denote the thermal average. The definition of the scalar order parameter (sometimes called nematisation) is deduced from the local second-rank order parameter tensor, M αβ (r) = S α r S β r − 1 2 δ αβ . After space averaging, the traceless tensor L −d r M αβ (r) admits two opposite eigenvalues ± 1 2 η corresponding to eigenvectors n + and n − . The order parameter density is defined after thermal averaging by M 2 (T ) = η . The associated susceptibility is defined by the fluctuations of the order parameter density, χ M 2 (T ) = 4L d k B T ( η 2 − η 2 ). Simulations are performed using a standard Wolff algorithm adapted to the expression of the nearest neighbour interaction [25,26]. The spins are located on the vertices of a simple square lattice of size L 2 with periodic boundary conditions in the two directions. Usually 10 6 equilibrium steps were used (measured as the number of flipped Wolff clusters) and 10 6 Monte Carlo steps for the evaluation of thermal averages (the autocorrelation time at k B T /J = 0.2, L = 16 is of order of 30 MCS, hence the numbers of MC steps corresponds roughly to 3.10 4 independent measurements for the smallest size). A first qualitative description of the behaviour of the system is provided by the temperature dependence of thermodynamic quantities [27]. The specific heat has a maximum which does not seem to increase substantially with the system size. This might be the sign of an essential singularity around a temperature k B T c /J ≃ 0.70. From the order parameter variation, a smooth transition is suspected, since there is no evolution toward a sharp jump. The susceptibility displays a non conventional behaviour at low temperature, increasing with the system size, which indicates a possible topological transition with a critical low temperature phase where the susceptibility diverges at any temperature.

Characterization of the low-temperature phase
We assume the existence of a critical phase at low temperatures as suggested by the temperature dependence results. The properties of the phase transition may be studied using i) Standard Finite-Size Scaling technique (FSS): in the critical low temperature phase of a model which displays a topological transition, the physical quantities behave like at criticality for a second-order phase transition, with power law behaviours of the system size. The difference is that in the critical phase, the critical exponents depend on the temperature and for any temperature below the transition one has e.g.
Here η M 2 (T ) denotes the correlation function critical exponent, ii) Rescaling of the density profiles (or "Finite Shape Scaling", FShS): conformally covariant density profiles or correlation functions are exepected at any temperature below the transition T c . They transform according to G(w 1 , w 2 ) = |w ′ (z 1 )| −xσ |w ′ (z 2 )| −xσ G(z 1 , z 2 ) through conformal mapping w(z) where w labels the lattice sites in the transformed geometry (the one where the computations are really performed), while z refers to the infinite plane where the two-point correlations take the standard power-law expression G(z 1 , z 2 ) ∼ |z 1 − z 2 | −ησ ), and x σ = 1 2 η σ is the scaling dimension associated to the scaling field under consideration. Rather than two-point correlation functions, it is even more convenient to work with density profiles m(w) in a finite system with symmetry breaking fields along some surfaces in order to induce a nonvanishing local order parameter in the bulk [28][29][30]. The density m(w) will be M 2 (r) = Q 2 (S r · h ∂Λ ) . In the case of a square lattice Λ of size L × L, with fixed boundary conditions along the four edges ∂Λ, one expects (details may be found e.g. in [29]) Another conformal mapping which has been applied to many two-dimensional critical systems is the logarithmic transformation w(z) = L 2π ln z = L 2π ln ρ+i Lϕ 2π which maps the infinite plane onto an infinitely long cylinder of perimeter L. The correlation functions along the axis of the cylinder (let say in terms of the variable u = L 2π ln ρ) decay exponentially at criticality, G(u 1 , u 2 ) ∼ exp[−(u 2 −u 1 )/ξ], the correlation length amplitude on the strip being universal, ξ = L πη [31]. Simulations at different temperatures in a system of size 10×10000 were performed and the correlation function exponent thus follows from the linear behaviour The resulting η−exponent deduced from these different techniques is shown in Fig. 3. The variation of η with the temperature is very similar to the one observed in the case of the XY model and confirms the QLRO nature of the LT phase of the model. The low temperature linear variation of this exponent is easily understood within the spin wave approximation. For O(2) model with nearest neighbour interactions described by arbitrary polynomial in S r · S r ′ , one is led to an effective harmonic Hamiltonian 1 2 J (r,r ′ ) l(θ r − θ r ′ ) 2 . It yields power-law correlations, with a decay exponent given by The comparison is made visible in the figure (with m = 2 and l = 128/11), and as expected, the lower the temperature the better the SW approximation.

Critical behaviour at the deconfining transition
Not only the low temperature behaviour of η is interesting, the precise value of the η exponent at the BKT transition where some deconfining mechanism should lead to the proliferation of unbinded topological defects is also of interest, since it really describes the universality class of the transition. An accurate value of the transition temperature is first needed. We performed a study of the crossing point of U 4 Binder cumulant for very large statistics (30 × 10 6 MCS) and large system sizes (squares of L = 64, 80, 96 and 128 with periodic boundary conditions). The results shown in Fig. 4 indicate a transition temperature of k B T BKT /J = 0.7226.
Then this temperature is used to perform Finite-Shape Scaling using the algebraic decay of density profiles inside a square with fixed bounday conditions. These simulations are time-consuming, since the autocorrelation time increases in the low temperature phase when T evolves towards the deconfining transition and a rather large number of Monte Carlo steps is needed to get a satisfying number of independent measurements. For sizes L = 16, 32, 64, we used 10 6 MCS for thermalization and 30 · 10 6 for measurements, while "only" 20 · 10 6 for the largest size 128. The exponential decay of two-point correlation functions along the torus cannot be applied at the BKT transition, since the system size being quite larger than in a square geometry, the number of MC iterations required is by far too large. In Fig. 4 we plot the "effective" exponent η eff (L) measured at T BKT for different system sizes as a function of the inverse size. An estimate of the thermodynamic limit value (L → ∞) can be made using a polynomial fit (the results of quadratic and cubic fits are respectively 0.118 and 0.122), but it is safer to keep the three largest sizes available, L = 32, 64, and 128, for which a linear dependence of η eff (L) with L −1 is observed. Taking into account the error bars, crossing the extreme straight lines leads to the following value for the correlation function exponent at the deconfining transition This value is essentially half the Kosterlitz value for the XY model.

Summary and open questions
The results obtained are essentially the following: -The P 4 O(2) model displays a BKT-like transition with QLRO in the LT phase where SWA nicely fits the nematisation temperature-dependent exponents η(T ) when T → 0.
-The value of the exponent at the deconfining transition is supposed to reach a universal value. The numerical estimate is close to 0.125.
We believe that our conclusions are safe in which concerns the existence of a QLRO phase at low temperature. The results may be understood through a naive comparison with clock model in 2d. Increasing the order of the interaction polynomial indeed increases the number of deep wells which stabilise the relative orientation of neighbouring spins, and one is thus led to a system which is quite similar to a planar clock model with a finite number of states, unless the fact that here we keep a continuous spin symmetry which prevents from any "magnetic" longrange order at finite temperature. The clock model is in the Potts universality class when q = 3, but at q ≥ 4, it displays a QLRO phase before conventional ordering at lower temperatures [12]. Combining this analogy with the requirements of Mermin-Wagner-Hohenberg theorem for continuous spin symmetry gives a natural explanation for our results. For any type of nearest-neighbour interaction (in P 1 , P 2 or P 4 ) the behaviour seems to be always described by a BKT-like transition. The similar observation that a two-component nematic model renormalises in two dimensions towards the XY model was already reported in Ref. [32]. The transition is likely driven by a mechanism of condensation of defects, like in the XY model, but due to the local Z 2 symmetry not only usual vortices carrying a charge 1 are stable, but also disclination points carrying charges 1/2 should be stable. The role of these defects might be studied in a similar way than in the recent work of Dutta and Roy [21], by the comparison of the the transition in the pure model and in a modified version where a chemical potential is artificially introduced in order to control the presence of defects. Now, the observed value of the exponent η at the transition temperature seems a bit strange. A naive application of the mechanism discussed in section 2 leads to a value two times larger. Apart from minor modifications, the same procedure applies. The Hamiltonian (15) is of the form (7) with K replaced by Kl which does not affect the universal properties, but changes the critical temperature by the same factor. The other modification comes from the local Z 2 symmetry in the nematic model. This changes the charges n from integers to half-integers, which modifies the factor of 2π in equation (11) to π/2. The divergence of the perturbative term in (12) then occurs at 1 2 πK c l = 4, where the spin wave exponent (22) becomes η c = m 2 /16 which takes the value 1 4 when m = 2. We thus have a huge discrepancy of 100 % between the numerical determination and this prediction! A probable source of this discrepancy may be in the mechanism involved. There are no negative charge in the nematic model, so the pair interaction has a different structure and a more refined analysis is therefore desirable.