57 | Study of phase transitions and structural order in perturbed nematic liquid crystals Študij faznega prehoda in strukturnega reda v perturbiranemu nematičnemu tekočemu kristalu Amid Ranjkesh * , Milan Ambrožič, Mitja Slavinec Faculty of Natural Sciences and Mathematics, University of Maribor, Koroška 160, 2000 Maribor, Slovenia E-mails: amid.ranjkesh@uni-mb.si, milan.ambrozic@uni-mb.si, mitja.slavinec@uni-mb.si * The corresponding author. Abstract: We studied the typical domain size and configuration character of randomly perturbed nematic liquid crystal system, or in general the system of rod-like objects which interact via a Lebwohl–Lasher-type interaction. We described their local direction with a headless unit director field. We introduced into the system the impurities of concentration p, which impose the random anisotropy field-type disorder to directors. We studied the domain-type pattern of molecules as a function of p, anchoring strength W between a neighboring director and impurity, temperature, and history of samples. In simulations, we quenched the directors either from the random or from homogeneous initial configuration. Our results show that a history of system strongly influences: i) the average domain coherence length; and ii) the range of orientational ordering in the system. In the random case, the obtained order is always short-range (SRO). On the contrary, in the homogeneous case, SRO is obtained only for strong enough anchoring W and large enough concentration p. In other cases, the ordering is either quasi long range (QLRO) or long range (LRO). We further studied field-induced memory effects for the random initial configuration. With increasing external ordering field, either QLRO or LRO is realized. This ordering is preserved even if the field is switched off, and its degree saturates for a large enough value of the field. Therefore, one can control the degree of global ordering and average domain coherence size by temporarily exposing the system to an external ordering field. Such systems could be exploited as soft matter based phase-memory devices. Key words: Lebwohl-Lasher model, nematic liquid crystals, disorder, memory effect, orientational order. Povzetek: V članku obravnavamo tipične konfiguracije naključno motenih nematskih tekočih kristalov in velikosti njihovih domen. Gre za sistem paličastih delcev, ki delujejo eden na drugega z interakcijo tipa Lebwoh-Lasher. Lokalno smer teh delcev opišemo z vektorskim (direktorskim) poljem. V sistem dodamo nečistoče s koncentracijo p, ki povzročajo naključno neurejenost direktorskega polja. Obravnavamo domensko strukturirano mrežo molekul kot funkcijo p, jakosti sklopitve W med sosednjima direktorjema in nečistočo, temperature in zgodovine vzorca. V numeričnih simulacijah smo začetno direktorsko polje dobili z nenadno ohladitvijo, ali pa s homogenimi začetnimi pogoji. Rezultati kažejo, da zgodovina sistema bistveno vpliva na: i) povprečno koherentno dolžino domene in ii) doseg orientacijskega urejanja sistema. V primeru naključne konfiguracije, dobimo vedno red kratkega dosega, v primeru homogenega urejanja pa dobimo red kratkega dosega le v primeru dovolj močne sklopitve W in dovolj velike koncentracije p. V ostalih primerih dobimo kvazi-dolgi doseg ali dolgi doseg interakcije. Obravnavali smo tudi spominske efekte naključne začetne konfiguracije na podlagi zunanjega polja. S povečevanjem zunanjega urejevalnega polja dobimo kvazi-dolgi doseg, ali pa dolgi doseg interakcije. Tako Izvirni znanstveni članek NARAVOSLOVJE ANALI PAZU 3/ 2013/ 2: 57-67 www.anali-pazu.si STUDY OF PHASE TRANSITIONS AND STRUCTURAL ORDER IN PERTURBED NEMATIC LIQUID CRYSTALS | 58 ureditev ohrani sistem tudi potem, ko zunanje polje odstranimo. To pomeni, da lahko kontroliramo stopnjo globalne urejenosti in velikost povprečne koherentne dolžine, če sistem začasno izpostavimo zunanjemu polju. Taki sistemi bi se lahko uporabili kot spominski sistemi na bazi mehkih snovi. Ključne besede: Model Lebwohl-Lasher, nematski tekoči kristali, nered, spominski efekt, orientacijska urejenost. 1. Introduction The phase behavior of fluid systems confined to porous hosts and of two-component mixtures is of wide interest in condensed matter physics, possibly because it stands at the crossroads of several issues of general importance, such as the effects of finite size and quenched disorder [1-6]. The isotropic to nematic phase transformation of thermotropic liquid crystals embedded in a number of porous matrices has been experimentally characterized using various techniques [7]. Two- component mixtures, such as nematic liquid crystals (NLC) with nano-particles (NP), can in general exhibit behaviors that are not encountered in either of isolated components, opening gates to new applications. Such systems are expected to play an important role in the emerging field of nanotechnology and in composites with extraordinary material properties. Great efforts have been devoted to the study of the thermodynamic, optical and dynamic properties of NLC merged in low-density silica aerogels [7]. NLC show relatively strong responses, even to local low-energy excitations, so their structure can be readily controlled by the confining surfaces and/or applying external magnetic or electric field. We confined our attention to thermotropic NLC in which liquid crystal phase is induced by lowering the temperature from the ordinary liquid (isotropic) phase. In the bulk nematic phase NLC molecules tend to be oriented homogeneously along a single symmetry breaking direction. At the mesoscopic level, the average local orientational ordering is commonly described by the nematic director fieldn  . The directions n   of this unit vector field are physically equivalent, reflecting the so called head-to tail invariance of NLC phase on the mesoscopic scale. For creating domain patterns, it is necessary to quench NLC from isotropic to the lower symmetry nematic phase [8]. A randomly chosen configuration of the symmetry breaking field n  is established in causally disconnected parts. This process is based on local fluctuation mediated preferences. Subsequently a domain structure appears which is well described by a single domain length  d. The reason behind this is the continuous symmetry breaking and causality. The basic features of domain pattern dynamics in a pure bulk are described by the Kibble- Zurek mechanism [9, 10]. It was originally introduced to explain the formation of topological defects in the early universe following the big bang [9]. Formation of topological defect at the domain wall is energetically costly due to high concentration of domain walls and defects. To reduce these costs domains grow with time what is enabled by mutual annihilation of defects [11, 12]. The homogeneous structure in the pure bulk system is gradually reached, but if impurities are present, they can pin the defects [13, 14]. Accordingly, the domain structure can be preserved. In this contribution, we use the Lebwohl-Lasher lattice model [15-17] for simulations of randomly perturbed NLC configurations. We show that in certain circumstances impurities can stabilize domain pattern giving rise to short range ordering in the nematic LC phase. 1.1. Model The three-dimensional (3D) spin model of NLC is represented by a rectangular simulation cell which is a lattice of M  N  L sites (typically, we choose M = N = L = 70, in order to make a compromise between the accuracy and calculation times). Each site is enumerated by a set of indices (i, j, k), where 1  i  M, 1  j  N and 1  k  L, and it is occupied by a spin S (i, j, k)  S ijk. The spin is a unit vector and may point in arbitrary direction (variant of the Heisenberg spin model). The neighboring NLC spins tend to align in parallel directions. In addition, we supposed that, instead of NLC, some randomly positioned sites represent impurities, which tend to reorient locally the spin orientation in random directions. The probability for the specific site to contain the impurity is p (which is shortly called impurity concentration) giving on average pMNL impurities in the lattice. Accordingly, the spin S represents either the NLC nematic director or the preferred direction of the impurity (i. e., the direction that the impurity tends to enforce to neighboring NLC spins), concerning whether the site contains NLC or impurity. The computer random number generator sets the static impurity positions and preferred ANALI PAZU, 3/ 2013/ 2, str. 57-67 Amid RANJKESH, Milan AMBROŽIČ, Mitja SLAVINEC 59 | directions in advance. In addition, a homogeneous field (electrical, magnetic, etc.) of size B in the arbitrary direction can also be applied to the cell that tends to align all the spins in its direction. The strengths of all three kinds of interactions are given by positive parameters J, W and B, which are called the spin coupling constant, the impurity anchoring strength and the field strength, respectively. The free energy functional is modeled by the sum of the terms of all spin sites:   ijk ijk f F , (1) Where the energy term f ijk consists of the “spin” and the “field” parts [18]: 2 2 2 , . . . . 1 ( ) ( ) 2 ijk ijk n n ijk n n ijk B i j i f J S S B S e       (2) The first part of the free energy functional includes the spin interactions between the six nearest neighbors (denotation n.n.). For convenience, the factor ½ is added in order not to count the interaction between each spin pair twice. The coupling constant J ijk,n.n. can have one of the three values: 1) 0, if both neighbors are impurities (since they are fixed and do not interact); 2) J, if both neighbors are NLC spins; and 3) W for mixed types. The second part of Eq. (2) means the interaction of the external field with the spins. The field B with fixed direction e B acts only on the sites with NLC. The free energy and the three interaction magnitudes (J, W and B) are next renormalized with respect to the coupling constant J by setting J =1. The equilibrium spin configuration is obtained by minimizing the total free energy with respect to all the spins. In order to satisfy the normalization of the spin vectors, S ijk 2 = 1, the total free energy must be rewritten as: * * ijk ijk Ff   (3) where: ijk ijk ijk ijk f S f    ) 1 ( * 2  , (4) The additional multiplication parameters  ijk are introduced, which also have to be found in order to solve the system. The total free energy (3) is minimized by setting to zero its derivatives with respect to 3(1  p)MNL nematic spin components (remember that impurity sites are intact). If all  ijk factors are first expressed (combining both minimization and spin normalization equations), we are left with the (1 p)MNL vector spin equations: 2 , . . . . .. ( , ) ( , ) 0 ijk ijk n n ijk n n ijk B nn R J g S S B g S e     , (5) Where the vector g is a function of two vectors, which is defined as:   2 1 2 1 2 1 2 1 ) ( ) ( ) , ( v v v v v v v v g              . (6) If the system equilibrium has not been reached yet, the left side R ijk of Eq. (5) is not zero and is called residuum. It is half the derivative of the free energy with respect to individual spin: ijk ijk S F R       2 1 . (7) If we are interested only in equilibrium of the system without thermal fluctuations (zero temperature), the system of equations (5) is solved with over-relaxation method where each NLC spin is corrected in proportion with the residuum R ijk. However, when thermal fluctuations of the spins have to be considered, the final equilibrium state of the system is calculated through the real-time relaxation process including thermal disturbances. The change of spin components in the time step t is: T ijk ijk ijk S old S F kT t D old S new S             ) ( ) ( ) ( . (8) Here D is the appropriate degenerated rotation diffusion tensor of the system and k is the Boltzmann constant. The first term on the right side of (8) corresponds to mechanical torque that tends to rotate the spin towards equilibrium, while the second term S T represents random thermal fluctuation. If we introduce dimensionless time step and temperature, t* = D t, T* = kT/J, and use Eq. (8) with dimensionless residuum vector, we obtain the corresponding numerical equation: T ijk ijk ijk S old R T t old S new S           ) ( * * 2 ) ( ) ( . (9) The appropriate dimensionless time interval t* is chosen as the input parameter according to literature: t*  0.016. In order to set the thermal fluctuation vector S T, it is convenient to consider first their rotation of the spins in their own frame (with the local z axis aligned with the spin direction). Only the rotation of the spin about perpendicular local x  and y axes is relevant and we suppose both rotations to be independent from each other. So we can write for the rotational extra energy: ) ( 2 2 2 1 y x J c E      , (10) Where  x and  y are the corresponding rotation angles and c 1 is the appropriate dimensionless constant (input parameter). In accordance with canonical STUDY OF PHASE TRANSITIONS AND STRUCTURAL ORDER IN PERTURBED NEMATIC LIQUID CRYSTALS | 60 distribution, we can write the infinitesimal probability for the range of angles within the phase-space element d  x  d  y: y x y x d d kT J c C dP                 2 ) ( exp 2 2 1 , (11) Where C is the normalization constant and k is the Boltzmann constant. Eq. (11) conveniently rewrote in the two-dimensional polar coordinate system:  d d T c C dP              * 2 exp 2 1 , (12) Here  is the magnitude of the vector (  x,  y) and the angle  defines its direction with respect to the local x- axis. To understand this, we may imagine the thermal fluctuation as the rotational kinetic free energy with two independent components of angular velocity  x and  y. Therefore, in the short time the spin is actually rotated by the angle  about the axis that is perpendicular to the local z axis and makes the angle  with the local x axis. While the distribution of  is uniform in the interval (0, 2 ), the probability distribution function for  in the interval (0, ) is:              2 2 2 exp ) (  C p , (13) with  2 = T*/c 1. The distribution function was simulated by our computer random generator. After the thermal rotation of the spin in its own frame the coordinates of the rotated spin must be transformed again to the laboratory system. The thermal rotation of the individual spin is finally written with the following vector transformation formula:                                            cos sin cos sin sin 0 1 1 1 1 1 2 2 2 2 2   z z y z x z z y x z y z z x T S S S S S S S S S S S S S S S  (14) In the matrix, there are components of the spin before rotation; it corresponds to two rotations of coordinate system; first for the spherical azimuthal angle about the z axis, second for the polar angle about the new y axis. This is not the only possibility for the transformation matrix, because it is only important that the laboratory z- axis rotated to match the spin direction (what can be seen from the third column of the matrix). The vector to the right of the matrix is the rotated spin in its own coordinate system. The orientational ordering of the spin system can be characterized by the traceless symmetric order parameter tensor with 3 3 components: I S S Q n ijk m ijk mn 2 1 2 3 , ,     , (15) where S ijk,m is the m-th component of the spin S ijk. The brackets <...> denote the average of the quantity through the simulation cell. I is the identity matrix. The scalar order parameter S is defined as the largest eigenvector of the matrix (15). The correlation function measures both short- and long-range spin orientation ordering. It is defined in a manner analogous to the order parameter tensor: 2 1 ) ( 2 3 ) ( 2      mnl ijk S S r G   , (16) where <...> is the statistical average of the squares of the scalar products of only those spin pairs which are separated by distance r. We have calculated it numerically after relaxation of the spin system. In the cases of SRO or QLRO it holds ( ) 0 Gr    . In the case of LRO it follows 2 () ~ G r g S     . In general, one expects an exponential decay towards a saturated value of G(r) on increasing r for both LRO and SRO. On the other hand, for QLRO algebraic decay of correlations () G r r    is expected [19]. In order to obtain structural details from G(r) for a finite system the correlation function is fitted using the empirical ansatz (0) 00 () kr G r a e b   (17a) or 1 (1) 1 1 () kr ae G r b r   (17b) where a 0, a 1, b 0, b 1, , k and k 1 are adjustable parameters. Note that it roughly holds 2 01 ~~ b b S , and 1/k   estimates an average linear size of a relatively well correlated region, referred to as a domain. The ansatz (0) () Gr or (1) () Gr is suitable for structures exhibiting either LRO or SRO. The parameters 2 02 ~~ b b S measure the LC ordering on large scale: S = 0 indicates SRO while S > 0 indicates LRO or QLRO. On the other hand, the expression (2) () Gr is more appropriate for studying structures possessing QLRO: (2) 2 2 () a G r b r   (17c) a  b  and   are adjustable parameters. One expects G   r → ∞) → 0. However, the decay of correlations with ANALI PAZU, 3/ 2013/ 2, str. 57-67 Amid RANJKESH, Milan AMBROŽIČ, Mitja SLAVINEC 61 | distance is relatively weak, and finite-size effects are expected to be important. The estimates for the power law coefficient  could be determined using a finite size scaling analysis [20] of S  SL   using the relationship : SL    (18) for SRO, the observed mean order parameter in a system of N particles would be given by S ∼ L -3/2 , implying  = 3 2 in this case. By contrast, in the case of LRO, the order parameter tends to a finite non-zero value for large system size L , and hence  = 0. QLRO is signalled by intermediate values of with  =  /2, so long as ≤ 3, although if  ≥ 3, the SRO finite-size scaling result is recovered. 2. Results and discussion Typical size of the lattice used in our calculations was M = N = L = 70. Before systematic variation of free parameters in our model (e. g., impurity concentration and anchoring strength), we adjusted the parameter c 1 in Eq. (10) so that the nematic scalar order parameter S drops to zero due to thermal fluctuations (for small concentrations of impurities) approximately at the value: kT/J  1. After several trials, we got roughly c 1 = 55 and we kept this value for all consequent simulations. Next, we chose two specific values of anchoring strength, W = 2 and W = 0.5, and varied systematically impurity concentration p and temperature T. For specific value of p, we started simulations with some low temperature (T = 0.1), set the desired initial spin configuration and then gradually increased temperature. Moreover, we tested the system behavior upon increasing temperature according to four different scenarios relating to initial system configuration: a) Random initial orientation (complete disorder) of nematic spins at the lowest temperature (T = 0.1) followed by spin’s relaxation. However, for each higher temperature, the numerical procedure takes the last (relaxed) spin configuration of the previous (a little lower) temperature before relaxation at this temperature begins. We shortly denote this scenario by R (random). b) Random initial orientation for all temperatures, not only the lowest one. We shortly denote this scenario by RR (random-random). c) Homogeneous (in one direction) initial alignment (complete order) of nematic spins at the lowest temperature (T = 0.1) followed by spins’ relaxation. However, for each higher temperature, the numerical procedure takes the last (relaxed) spin configuration of the previous (a little lower) temperature before relaxation at this temperature begins. We shortly denote this scenario by H (homogeneous). d) Homogeneous initial orientation for all temperatures, not only the lowest one. We shortly denote this scenario by HH (homogeneous- homogeneous). We also note that in any case, the configuration of impurities is unchanged during rising the temperature. We also took care that for some given value of p, the configuration of impurities was the same for the four scenarios described above. Experimentally scenario H means that in the beginning we used very strong homogeneous external field only at the lowest temperature, to align all spins. HH means that the initial spin configuration is aligned with the strong field for every temperature before spin relaxation. R means that the system is frozen suddenly from isotropic phase to the lowest temperature, and then the temperature is slowly increased. RR means the system is suddenly frozen from isotropic phase for every considered temperature before spin relaxation. In simulation, there was no external field (except for considering very strong initial field to achieve aligned initial configuration, which is then turned off). We followed the evolution of the scalar order parameter S and  belonging to the correlation function (16) according to the form (17a). Since both impurities and thermal fluctuations contribute to orientational disorder in the nematic phase, both parameters S and  are diminished when either impurity concentration or temperature are increased. For each temperature, we made 5000 sweeps (cycles) over the entire lattice. This was in most cases enough to obtain qualitatively correct temperature dependences of the macroscopic parameters. The large numbers of sweeps would require much more computational time. For comparison with another research group [21] used 10 5 cycles for their analysis of memory effects in nematics with quenched disorder although they commented that this number was far larger than necessary to retain high precision. In order to eliminate small thermal fluctuations of the macroscopic parameters for each temperature, we averaged them for the last 40 sweeps out of 5000. Nevertheless, figures still have evident small statistical irregularities in the temperature dependences of calculated parameters. In order to check quantitatively how well the 5000 sweeps present the thermal equilibrium structure, we followed the change of some macroscopic variables with increasing number of sweeps. STUDY OF PHASE TRANSITIONS AND STRUCTURAL ORDER IN PERTURBED NEMATIC LIQUID CRYSTALS | 62 We did this for a few different sets of parameters (temperature, impurity strength and concentration, and initial condition of spins – random or aligned). Typical situation is presented in Figures 1 3. The parameters for all these figures are W = 2, p = 0.2, T = 0.1 and HH scenario. Figure 1. Dependence of the free energy F (calculated per site) on the number of sweeps. Figure 2. Dependence of the order parameter S on the number of sweeps. Figure 3. Dependence of the coherence length  on the number of sweeps. 2.1. Different Histories In Figure 4 we plotted typical obtained G(r) profiles, with r normalized to nearest-neighbor distance a 0. This is the simplest thing which can be done now since a 0 is already used in Figure 4. In the case of short-range order G(r) drops to zero for 𝑟 / 𝑎 0 >> 1. On the other hand, for QLRO or LRO the correlation function reaches a finite value in the limit 𝑟 / 𝑎 0 >> 1. Figure 4. G( 𝑟 ) profiles. ○: random history, ×: homogeneous history. 𝑁 = 70, 𝑊 = 1, 𝑇 = 0.5, p = 0.1. ANALI PAZU, 3/ 2013/ 2, str. 57-67 Amid RANJKESH, Milan AMBROŽIČ, Mitja SLAVINEC 63 | In Figure 5 we considered the value of S as a function of temperature T, we used the scenarios H and HH since this parameter for random initial spin configurations (R and RR) are zero. The values of this parameter for HH are higher than for H scenario as expected. The dependence of the correlation length  on temperature T, spin concentration and scenario is plotted in Figure 6 for comparison. It is evident from Fig. 6 that the correlation length  is not significantly different for the four scenarios, except for smaller values of impurity concentrations p. However, it can be noted that correlations length for R is systematically higher than for RR, while that for H it is systematically higher than that for HH. Rearranging spins (either aligning or randomizing) for each temperature separately obviously lowers the correlation length a little. We showed only the situation for W = 2 since the lines for W = 0.5 are much too irregular (large fluctuations) for the picture to be representative. It is evident from Fig. 6 already that the fluctuations in  are larger for smaller impurity concentrations. Averaging over several repetitions of the same temperature scan would smooth these curves but this requires much more computational time. Figure 5. The degree of orientational ordering as a function of T for different values of p (2.5 % to 20 %); (a) W = 2 and (b) W=0.5. Scenarios: H (full forms) and HH (empty forms). Figure 6. Dependence of correlation length  as a function of T in different P (5 % to 20 %) and, W = 2. Scenarios: H (●), HH (○), R ( ■) and RR (□). We made a brief study of temperature hysteresis behavior in the presence of impurities p (without external field). We took an example with W =1 and different value of p. We started with the aligned spin configuration at the starting temperature, and then the temperature was STUDY OF PHASE TRANSITIONS AND STRUCTURAL ORDER IN PERTURBED NEMATIC LIQUID CRYSTALS | 64 gradually decreased and next increased again. The results are shown in figure 7. Figure 7. History-dependence of the random system. Order parameter temperature-dependence S(T). “Up” (empty symbols) corresponds to increasing, “Down” (full symbols) corresponds to decreasing W = 1, N = 70. We next consider the range of orientational order within samples. For this purpose, we started from two significantly different histories of systems. We either initiated the simulations from LC structures homogeneously aligned along a single symmetry- breaking direction or from randomly aligned structures. We henceforth refer to these initial configurations as homogeneous and random histories, respectively. We first calculated structures originating from the random history. For all the studied cases, we obtained SRO. The G(r) profiles for p = 0.3 are shown in Figure 8. It can be concluded that with increasing temperatures T the domain size decreasing and SRO is observed. Figure 8. G(r) dependence for random histories. p = 0.3, 𝑊 = 1, 𝑁 = 70. In addition, we have made some finite size analysis of the system. We tested the resulting structure ordering for random initial configuration and for growing size of the system (M = N = L): from 20 to 80. For each N, we made 8 repetitions of the calculations (for the same system parameters), because of randomness of the impurities positions and orientations results in the distribution of all macroscopic parameters. Therefore, the results represented in Figure 9 are the averages of 8 repeated calculations. We also chose two representative temperatures: kT/J = 0.1, kT/J = 0.5 and two different boundary conditions, so called free and periodic. It seems for small size N, one suspects that the effects of the finite system become relevant. More specifically, the periodic boundary should result in larger value of order parameter S since the opposite faces of the simulation cube are connected. The analysis confirmed an expectation what can be seen on Figure 9. The Fitting of logarithmic graphs to linear functions showed the following values of the exponents in the equation (18) dependences:  = 1.77 for kT/J = 0.1 and free boundary,  = 1.85 for kT/J = 0.1 and periodic boundary,  = 1.88 for kT/J = 0.5 and free boundary,  = 2.13 for kT/J = 0.5 and periodic boundary. It seems that our simulations correspond more to SRO than to QLRO. Figure 9. Dependence of scalar order parameter S on the system size M = N = L; other parameters: p = 0.1,W = 2, kT/J = 0.1 and 0.5; RR scenario. ANALI PAZU, 3/ 2013/ 2, str. 57-67 Amid RANJKESH, Milan AMBROŽIČ, Mitja SLAVINEC 65 | Figure 10. Dependence of correlation length  on the system size M = N = L; other parameters: p = 0.1, W = 2, kT/J = 0.1 and 0.5; RR scenario. Figure 10 presents the dependence of correlation length  on system size. For both free boundary and periodic boundary condition (N) is slightly increasing function of system size, and is larger for periodic boundary. 2.2. Memory Effect Finally, we checked the memory effects of the system in connection with homogeneous external field, e. g. magnetic field. We compared the residual scalar order parameter S (after switching the field off) in two different procedures: zero-field-cooling (ZFC) and field cooling (FC). In both cases, we started with high temperature (above nematic-isotropic transition) and gradually lowered the temperature to some value deep in the nematic phase. We chose W = 1 and p = 0.14 for static impurities. The initial and the final temperature was T i = 1 and T f = 0.2, respectively, with the temperature step T = 0.05. In ZFC procedure the temperature was first gradually lowered from T i to T f with B = 0 (allowing the system to relax at each intermediate temperature). Then at T f = 0.2 the field was suddenly raised (in one step) from zero to some specific value B and the system was relaxed at B. Finally, B was suddenly switched off and the system relaxed at zero field. In FC procedure the field was suddenly raised (in one-step) from zero to B already at the temperature T i = 1 and the system was relaxed at B. Then keeping a field the temperature was lowered from T i to T f (allowing the system to relax at each intermediate temperature). Finally at T f = 0.2 the field was suddenly switched off, and the system relaxed at zero field. For both procedures, we started with random initial configuration and then followed the R scenario (the system remembered intermediate configurations); we also stored initial random configuration on file, so that it was identical for ZFC and FC case. The number of cycles was 10 4 for the single field step from 0 to B or opposite, while we took only 2000 cycles for one temperature step T to spare computational time (but since there are 16 temperature steps of 0.05 from T i to T f, this makes 32000 steps between the two temperatures). The size of the system was M = N = L = 70. We varied the field B to see its influence on the difference of the order parameter for both procedures. The results are presented in Figure 11. Figure 11. Dependence of residual S after removing the field on the field value for ZFC and FC procedures; p = 14 %, W = 1, T f = 0.2, R scenario. ●: ZFC, ○: FC. Briefly, in the study of memory effect, for some relatively small value of the field magnitude, after switching the field off, the FC procedure results are in higher value of residual order parameter than the ZFC procedure. 3. Conclusions We studied rod-like objects, namely nematic liquid crystals, within a cubic lattice, interacting via a Lebwohl– Lasher-type interaction. The structure of the system was STUDY OF PHASE TRANSITIONS AND STRUCTURAL ORDER IN PERTURBED NEMATIC LIQUID CRYSTALS | 66 described in terms of the director field S i, where the unit vector S exhibits head-to-tail invariance. Moreover, we introduced impurities of concentration p, which impose the so-called random anisotropy field. We analyzed the domain-type pattern of molecules as a function of p, anchoring strength W between a neighboring director and impurity, as well as initial configuration or history of the system. In random history, we quenched randomly distributed nematic spin orientations to a low temperature T << T  , and then allowed the system to equilibrate. In homogenous history, by contrast, we let a set of initially aligned nematic spins to relax. Our simulations showed that in both cases, for homogeneous and random histories, the order parameter S is reduced with increasing anchoring strength W and impurity concentration p, as well with increasing temperature. Similar conclusion holds for the correlation length. Finite-size scaling methods [20, 22] have been used to distinguish phases between LRO, QLRO and SRO. Our study regarding to finite-size scaling behavior in the N-I transition temperature for random history in different temperatures and two different boundary conditions reveals that in all cases the orientational order is SRO than QLRO. But for low values of p and W we can see QLRO. On the other hand, for homogenous history SRO appears only for very high value of anchoring strength W and impurity concentration p. We also analyzed field-induced memory effects for saturated configurations in the random initial configuration. Our study divulges that even a small external homogeneous field sets up LRO In the system. Acknowledgement A. Ranjkesh Acknowledges the Slovenian Research Agency (ARRS) for training and financing of young researchers. References 1. Zurek, W. H. Cosmological experiments in superfluid helium? Nature 1985, 317, 505. 2. Bray , A. J. Theory of phase-ordering kinetics. Adv. Phys 1994, 43, 357. 3. Kibble , T. W. B. Topology of cosmic domains and strings . J. Phys. A 1976, 9, 1387. 4. Imry, Y; Ma , S. Random-Field Instability of the Ordered State of Continuous Symmetry. Phys. Rev. Lett 1975, 35, 1399. 5. Feldman , D. E. Quasi-long-range order in nematics confined in random porous media. Phys. Rev. Let 2000, 85, 4886. 6. Repnik, R; Ranjkesh, A; Simonka, V; Ambrozic, M; Bradac, Z; Kralj, S. Symmetry breaking in nematic liquid crystals: analogy with cosmology and magnetism .J. Phys: Condens. Matter 2013, 25, 404201. 7. Crawford, G. P. ; Žumer, S. Liquid Crystals in Complex Geometries Formed by Polymer and Porous Network; Oxford University Press: London, 1996. 8. Bellini T.; Buscaglia, M.; Chiccoli, C. Nematics with quenched disorder: What is left when long range order is disrupted? Physical Review Letters 2000, 85(5), 1008. 9. Radzihovsky, L.; Toner , J. Anomalous Elasticity of Disordered Smectics.Phys. Rev. Let 1997, 79, 4214. 10. Popa-Nita,V ; Statics and Kinetics at the Nematic- Isotropic Interface in Porous Media, Eur. Phys. J 1999, 83, 12. 11. Popa-Nita, V.; Romano, S.; Nematic-Smectic A Phase Transition in Porous Media. Chem. Phys 2001, 91, 264. 12. De Gennes, P.G; Prost, J. The Physics of LiquidCrystals ; Publisher: Oxford University Press: Oxford, UK, 1993. 13. Virga, E.G . Variational Theories for Liquid Crystals. Publisher:Chapman Hall: London, UK, 1994. 14. Ambrožič, M.; Kralj, S.; Virga, E.G . Defect- enhanced nematic surface order reconstruction , Phys. Rev. E 2007, 75, 031708 . 15. Lebwohl, P.A; Lasher, G. Nematic-Liquid-Crystal Order-A Monte Carlo Calculation, Phys. Rev. A 1972, 6, 42. 16. Krasna, M; Cvetko, M; Ambrozic, M. Symmetry breaking and structure of a mixture of nematic liquid crystals and anisotropic nanoparticles, Beilstein J. Org. Chem 2010 , 6, No. 74. 17. Ranjkesh, A.; Ambrozic, M. ; Cordoyiannis, G.; Kutnjak, Z.; Kralj, S. History-Dependent Patterns in Randomly Perturbed Nematic Liquid Crystals . Advances in Condensed Matter Phys 2013, 2013, 505219. 18. Romano, S.; Computer simulation study of a Nematogenic Lattice-Gas model with fourth- rank interactions, Int. J. Mod. Phys. B 2002, 16, 2901. ANALI PAZU, 3/ 2013/ 2, str. 57-67 Amid RANJKESH, Milan AMBROŽIČ, Mitja SLAVINEC 67 | 19. Kutnjak, Z.; Kralj, S.; Lahajnar, Žumer, G. S. ; Calorimetric study of octylcyanobiphenyl liquid crystal confined to a controlled-pore glass. Phys. Rev. E 2003, 68, 021705. 20. Eppenga, R.; Frenkel, D. Monte Carlo study of the isotropic and nematic phases of infinitely thin hard platelets. Mol. Phys 1984, 52, 1303. 21. Serra F.; Vishnubhatla, K.; Buscaglia, M.; Cerbino, R.; Osellame, R.; Cerullo . G .; Bellini, T. Topological defects of nematic liquid crystals confined in porous networks . Soft Matter 2011, 7, 10945. 22. Bradac Z.; Kralj, S.; Zumer, S. Molecular dynamics study of the isotropic-nematic quench. Phys. Rev. E 2002, 65, 021705.