History-Dependent Patterns in Randomly Perturbed Nematic Liquid Crystals

We study the characteristics of nematic structures in a randomly perturbed nematic liquid crystal (LC) phase. We focus on the impact of the samples history on the universal behavior. The obtained results are of interest for every randomly perturbed system exhibiting a continuous symmetry-breaking phase transition. A semimicroscopic lattice simulation is used where the LCmolecules are treated as cylindrically symmetric, rod-like objects interacting via a Lebwohl-Lasher (LL) interaction. Pure LC systems exhibit a first order phase transition into the orientationally orderednematic phase atT = Tc on lowering the temperatureT.Theorientational ordering of LC molecules is perturbed by the quenched, randomly distributed rod-like impurities of concentration p. Their orientation is randomly distributed, and they are coupled with the LC molecules via an LL-type interaction. Only concentrations below the percolation threshold are considered. The key macroscopic characteristics of perturbed LC structures in the symmetrybroken nematic phase are analyzed for two qualitatively different histories at T ≪ T c . We demonstrate that, for a weak enough interaction among the LCmolecules and impurities, qualitatively different history-dependent states could be obtained.These states could exhibit either short-range, quasi-long-range, or even long-range order.


Introduction
Domains often appear in phases of broken continuous symmetry [1,2].The reason behind this is causality that is a finite speed of the information propagation [3].Consequently, in a fast enough phase transition in a different part of systems a different gauge component of the order parameter field is selected.This results in the appearance of regions in which the gauge field is essentially spatially homogeneous.These regions are referred to as domains and the mechanism of their formation as the Kibble-Zurek mechanism [4].At boundaries between domains topological defects (TDs) [5] reside carrying a topological charge .The total topological charge  tot of the system is conserved.At a given time, a domain pattern is well determined by a single characteristic size   .The initial size of visible domains (the so-called protodomains) is dictated by the phase transition quench rate [4,6,7].The domain walls are energetically costly, and, consequently, the average domain size grows with time  following the scaling law   ∝   , where  is a universal scaling coefficient [1].This growth is enabled by the annihilation of defects and antidefects.In ideal pure systems, a single domain would be gradually formed as a function of time.However, when impurities are present, they could stabilize the domain pattern depending on their concentration  and interaction with the host [2].In some cases, impurities could be the carriers of topological charge, and thus, they give rise to additional TDs in order that  tot is conserved.The resulting domain pattern might be very complex and exhibit strong memory effects [8].Different domain configurations could lead to effectively new materials [9], which can be exploited in various electrooptic applications, memory elements, and so forth.Consequently, it is of interest to understand in detail the mechanisms that control the key domain pattern characteristics.
Nematic liquid crystals (LCs) comprise a convenient, simple, and experimentally relatively easily accessible system to study the domain patterns.Due to their soft and liquid-like Advances in Condensed Matter Physics character, one can tune their structure using a relatively low power.In addition, they are optically transparent and anisotropic, which allows one to monitor the LC patterns using simple optical experiments (e.g., observing the textures by means of polarizing optical microscopy).For simplicity we henceforth consider LCs formed by rod-like molecules.In thermotropic LCs, a nematic () phase is reached on lowering temperature from the isotropic () phase.The latter phase is characterized by an isotropic liquid character.In the nematic phase the continuous symmetry is broken, and the molecules tend to be on average aligned parallel along a common direction, referred to as the director field ⃗ , where the states ± ⃗  are equivalent and | ⃗ | = 1.In equilibrium bulk  phase ⃗  is homogeneously aligned along a single symmetrybreaking direction.
In the past few decades, several studies have analyzed the phase transition and structural behavior of mixtures of LCs with aerosil spherical particles [10][11][12][13][14].These mixtures represent adequate experimental model systems to study the influence of qualitatively different origins of disorder on phases with broken continuous symmetries.By changing the mass density of the aerosil particles, at least three qualitatively different regimes are established [11,12].Various studies report that either short-range order, quasi-long-range order or even long-range order might be established.In some cases glassy-like behavior and pronounced memory effects are also reported [10][11][12][13][14].
In this paper, we focus on the impact of history on nematic structures in randomly perturbed LCs.We use random anisotropy nematic (RAN) lattice model in which rod-like impurities and LC molecules interact via a Lebwohl-Lasher type interaction [15].Such systems roughly simulate the mixtures of LCs and aerosils.The macroscopic nematic structural properties are studied as a function of the impurities concentration , for a given interaction strength  between impurities and LC molecules, and for two extremely different sample histories.It is shown that the nematic structures can exhibit either long-range order (LRO), quasi-longrange order (QLRO), or short-range (SRO) order on varying the disorder strength and sample history.In the case of QLRO orientational correlations decay algebraically with distance.Topological reasons for random-field type of disorder in mixtures of LCs and impurities are described in Section 2. The model we use to study macroscopic properties of systems of our interest is defined in Section 3. The corresponding results are presented in Section 4. Key conclusions are summarized in the last section.Some supporting derivations are assembled in the appendices.

Topology and Random-Field-Type Disorder
In mixtures of LCs with different particles (colloids or nanoparticles), to which we refer as impurities, the latter often introduces into the system a random-field (RF) type of disorder.In this subsection we present the basic mechanisms that mediate such a behavior.
At the continuum level, a local average orientation of LC molecules in the nematic phase is given by the director field ⃗  and the scalar order parameter  describing the amount of fluctuations about that direction.In the isotropic phase it holds  = 0, whereas a perfect alignment corresponds to  = 1.The equilibrium orientational configuration at a given temperature, for a fixed volume of the system, minimizes the free energy [16] The condensation free energy contribution term   determines the degree of ordering in the bulk LC in the absence of any elastic distortions.The elastic term   determines the elastic penalties if the system exhibits spatial inhomogeneities.The field term   describes the impact of the external electric or magnetic fields.The interfacial contribution presents the interactions at the LC-particle interfaces.Within the lowest order approximation these terms are typically expressed as [16] Here ,  * , ,  represent material constants,  and  are the representative elastic constants, ⃗  stands for an external ordering electric or magnetic field, Δ is the corresponding field anisotropy constant, ⃗ ] is the local normal of an LCimpurity interface, and  is the surface interaction constant.It roughly holds that the quantities , , , , Δ, and  are temperature independent and  ∼  2 [16].
In the bulk nematic equilibrium and in the absence of impurities, the LC orientational ordering is homogeneous along a single breaking direction.However, due to the broken continuous symmetry the nematic director field could exhibit topological defects (TDs), at which ⃗  is not uniquely defined.Their key feature is defined by the topological charge  [5,17,18].
The topological charge plays the key role in the classification of TDs in condensed matter.For a general -dimensional unit vector field ⃗  () = ( where  1 ,  2 , . . .,  −1 are the coordinates determining the ( − 1) dimensional sphere enclosing the defect and Π () is the normalization constant.For example, for  = 2 it holds ⃗  (2) = ( 1 ,  2 ), Π (2) = 2, and integration is performed along a closed line: For  = 3 it holds ⃗  (3) = ( 1 ,  2 ,  3 ), Π () = 2, and integration is performed along a closed surface: The stability of an isolated defect is guaranteed by the conservation of its topological charge.For example, a defect exhibiting a radial distribution of ⃗  () corresponds to  () = 1.Furthermore, the total topological charge of the vector field enclosed within the ( − 1) dimensional surface is conserved.The conservation laws of TDs are analogous to the laws of electric charges.They regulate the decay, merging, annihilation, and/or mutual transformations of TDs.
We henceforth focus our interest on the nematic LC ordering.Note that the nematic LC phase is characterized by the head-to-tail invariance of the director field ⃗ .Therefore, for  = 3 the sign of  does not play a role.For later convenience we assign to the nematic defects (antidefects) a positive (negative) value of .One typically refers to the nematic TDs bearing  = 1 ( = −1) as the nematic monopole or hedgehog (antimonopole, antihedgehog).The nematic monopole and antimonopole attract each other and they tend to be mutually annihilated into a defect-free state characterized by  = 0.The total topological charge  tot within a given volume is conserved if the boundary conditions remain unaltered.For example, for a homeotropic alignment at an enclosing surface (i.e., ⃗  is aligned along the local surface normal) it holds  tot = /2.The quantity  determines the Euler characteristic [5,17] of the closed surface and is topologically invariant (it does not change under smooth deformations).For example, for the spherical topology it holds  = 2 and  tot = 1.
In order to illustrate the impact of topology on the structural characteristics of the systems of interest, we consider a mixture of nematic LC and spherical impurities of radius .We set that the LC-impurity interfaces enforce a homeotropic anchoring (i.e.,  < 0 in (2d)).In the case of weak anchoring realized for /() < 1, the elastic term prevails over the interface tendencies.Consequently, the impurities have mainly a dilution impact that reduces the effective LC-LC interactions as illustrated in Figure 1(a).In case that the system boundary conditions do not enforce a topological charge, the total topological charge equals to  tot = 0.In the strong anchoring limit, realized for /() ≫ 1, each impurity effectively acts as a radial hedgehog bearing a topological charge  = 1.The topological charge conservation law requires that a compensating antidefect with  = −1 is created in the LC matrix.For appropriate LC elastic properties the combination of defectanti-defect forms an effective topological dipole [19,20] as shown in Figure 1 Figure 1: (a) In weak surface interaction limit an impurity does not essentially disturb LC environment.(b) Schematic presentation of topological defect dipole.The dipole is enforced by a spherical particle enforcing homeotropic anchoring at the LC-particle interface.Consequently, from LC perspective the particle carries topological charge  = 1 and can be viewed as a virtual topological defect.Due to the conservation of the total topological charge, an antidefect appears in the vicinity.The defect and antidefect attract each other but cannot mutually annihilate due to the virtual character of the defect.Note that the particle in absence of LC does not enforce any local orientational preference.However, topologically enforced formation of anti-defect necessarily breaks the initial isotropic symmetry.
a local orientational anisotropy.On the contrary, an isolated impurity does not enforce any anisotropy.In some cases the interactions between impurities (e.g., in mixtures of LCs and aerosils where the latter tend, to self-assemble; see [11,12]) give rise to an essentially randomly shaped network of impurities.Consequently, the orientational distribution of topological dipoles could be essentially randomly distributed within a system, imposing to the surrounding LC medium Advances in Condensed Matter Physics a kind of random field.These are the cases that are considered henceforth.

Theoretical Background
We use a semimicroscopic lattice model consisting of rodlike LC molecules and rod-like impurities.The orientational ordering of either of them at an th lattice site is presented by a unit vector ⃗   exhibiting head-to-tail invariance (i.e., states ⃗   and − ⃗   are equivalent).The LC molecules ( ⃗   = ⃗   ) and impurities ( ⃗   = ⃗   ) are henceforth referred to as LC spins and impurity spins, respectively.They occupy  3 lattice sites of a 3D cubic simulation cell.A typical value of  = 60 is chosen.An th site is occupied either by an LC or an impurity spin.Impurities of concentration  are randomly distributed; thus, the system possesses  3 impurities.Their orientations are also randomly distributed and quenched with time.

Interaction Potential.
The total interaction of the system is given as the sum over the lattice site contributions  = ∑    , where [13,[21][22][23]] The 1st term describes the short-range interaction   between neighboring spins.The coupling constant can have one of the three values: (1) 0, if both neighbors are impurities (i.e., their orientations are quenched); (2)   =  > 0, if both neighbors are LC spins; and (3)   =  > 0 for neighboring LC and impurity spins.The second part of (6) describes the interaction between the external field ⃗  =  ⃗   and the LC spins, where | ⃗   | = 1.At the mesoscopic level the first term in (6) is described by (2a), (2b), and (2d).Furthermore, the 2nd term corresponds to (2c).
In simulations we take into account the normalization | ⃗   | = 1 by introducing the Lagrange multipliers   in the functional  * = ∑   *  , where The functional  * is minimized with respect to LC spins ⃗   , yielding (1 − ) 3 equations: Here  ⇀   is commonly referred to as the residuum, and In the absence of thermal fluctuations (i.e., at  = 0), (4) is solved using the standard overrelaxation method.For finite temperatures the final (equilibrium or metastable) states are reached via the "real-time" relaxation process.The change of spin components in a time step Δ is given by Here  is the appropriate degenerated rotation diffusion constant of the system, and   is the Boltzmann constant.The second term on the right hand side of (10) corresponds to the mechanical torque which tends to rotate an LC spin towards equilibrium, while the third term Δ  ⇀   represents the random thermal fluctuations.Its determination is described in Appendix A.
We henceforth set  = 1 and introduce a dimensionless time step Δ * = Δ and temperature  * =   /.In the simulations it is set Δ * ≈ 0.01 in order to obtain sensible reference results in bulk nematic phase.
It follows Using this method the fixed point configurations are calculated, for which the macroscopic properties of the system do not anymore change as a function of time.Therefore, these configurations correspond either to equilibrium or to metastable states surrounded by relatively high energy barriers.

Measured Quantities.
From reached fixed point configurations the average tensor order parameter  of a system and the orientational correlation function () are calculated.The traceless symmetric order parameter tensor, describing the average behavior of the whole system, is defined as The brackets ⟨ ⋅ ⋅ ⋅ ⟩ denote a spatial averaging,  is the identity matrix, and ⊗ stands for the tensor product.The average scalar order parameter of the system  is defined as the largest eigenvector of .The correlation function is calculated as where ⟨ ⋅ ⋅ ⋅ ⟩ is the average over only those LC spin pairs which are separated by the distance .
In cases of SRO or QLRO it holds ( → ∞) → 0. In case of LRO it follows ( → ∞) →  ∞ ∼  2 .In general, one expects an exponential decay towards a saturated value of () on increasing  for both LRO and SRO.On the other hand, for QLRO algebraic decay of correlations () ∝  − is expected [21].
In order to obtain structural details from () for a finite system the correlation function is fitted using the empirical ansatz or where  1 ,   estimates an average linear size of a relatively well-correlated region, referred to as a domain.The ansatz  (1) () is suitable for structures exhibiting either LRO or SRO.On the other hand, the expression  (2) () is more appropriate for studying structures possessing QLRO.

Simulation Result
We first consider the typical temperature behavior of systems of interest using our modeling.Note that similar studies have already been performed using different approaches.Therefore, we just summarize the main results which are in line with these previous findings.The temperature behavior () of the degree of orientational ordering in presence of impurities, upon increasing and decreasing , is shown in Figures 2 and 3, respectively.In Figure 2 the () dependence is shown for concentrations  = 0.025,  = 0.05, and  = 0.075.One sees a hysteretic behavior related to the 1st order character of the transition.Therefore, the coupling strength between LC and impurity spins is below the critical value.Note that, due to the local ordering tendency of impurities, the isotropic phase is replaced by a paranematic phase for  > 0. In the presented cases the interaction between impurity and LC spins is weak enough so that the phase transition is preserved.In Figures 3(a compare the () variations for different values of  on increasing and decreasing , respectively.With increasing  the degree of ordering gradually decreases due to the increasing role of randomness.Simulations reveal that finitesize effects are still evident despite using relatively large samples.This is evidently shown in the plot  * * () in Figure 4, where  * * roughly estimates the temperature separating the paranematic and nematic phase on increasing .Impact of  on  * * for different interaction strengths  is shown in Figure 5.One sees a roughly linear  * * dependence on , where  * * decreases on increasing .Note that simple dilution mechanism predicts a linear dependence as it is demonstrated in Appendix B. Therefore, in the above-shown cases the impact of disorder does not qualitatively change dilution-driven  * * () behavior.We next consider the range of orientational order within samples.For this purpose we start from two significantly different histories of systems.We either initiate 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.In Figure 6 we plot typical obtained () profiles.In the case of short-range order () 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.
We first calculate structures originating from the random history.For all the studied cases we obtain SRO.The () profiles for  = 0.1 and  = 0.3 are shown in Figures 7(a) and 7(b), respectively.The size of domains  is estimated using the ansatz of Equation (14a), where  ∼ 1/ and is plotted in Figure 8.As expected with increasing value of  the size of domains shrinks.Note that with increasing the value of  one observes increasingly pronounced () dependence.With decreasing  the domains appear to be larger due to the increasing importance of elastic forces that tend to establish a spatially homogeneous structure.The observed SRO is in line with the Imry-Ma theorem [24], claiming that even infinitesimally weak RF-type disorder can destroy the LRO of a pure system reached via a continuous symmetry-breaking transition.
Next we consider structures obtained from the homogeneous histories.In this case, either LRO, QLRO, or SRO is obtained depending on values of  and .For low enough values of  and  one obtains LRO.In this case ( ≫  0 ) saturates at a finite value independent of , as shown in Figure 9. On increasing  or , quasi-long-range order is formed.In this case () dependence becomes pronounced as it is demonstrated in Figure 10.The linear dependence in a log-log diagram reveals a power law () dependence.The key features of the established QLRO are obtained by using the ansatz Equation (14b).In the weak coupling limit, LRO is formed with relatively weak variations in ().Consequently,   ∼ 0 is expected.On the other hand, in the strong coupling limit one expects in the thermodynamic limit that () = 0, that is, ( ≫  0 ) ∼  2 2 ∼ 1/ 3 for a finite number .Taking into account / 0 ∼ , where  corresponds to the linear simulation cell size, one expects  ∼ 3 in the strong interaction limit.
In Figure 11  is plotted as a function of  for  = 1 deep in the nematic phase.In accordance with our expectations we obtain  ∼ 0 for  ∼ 0, whereas  roughly approaches towards 3 on increasing  above the percolation threshold.The temperature dependence of  is analyzed in Figure 12.One sees that  increases with increasing .Such a behavior is expected because the thermal fluctuations are increasing, and thus, the same happens with the degree of disorder.From the diagrams we infer that  is a nonuniversal quantity, because it exhibits a dependence on  and material properties of the systems.For even higher values of coupling SRO is established.This is demonstrated in Figure 13 where  is plotted as a function of  for different values of .One sees that for  = 2 and  > 0.5 the disorder is strong enough to form SRO.

Conclusions
Our study is motivated by studies in magnetism which have shown that random-field-type disorder might give rise to different orientational structures [24,25].Because the main mechanism behind observed is continuous symmetry breaking, similar features are expected to be observed in liquid crystals.For this purpose we carried out systematic study of the impact of the samples history on the degree of orientational ordering in randomly perturbed nematic LCs.The dynamic Lebwohl-Lasher-type lattice model has been used.The systems of interest simulate to good extent the mixtures of LCs and aerosils or mixtures of similar type.In order to demonstrate the possible structural variation range that one could encounter, two limiting histories of samples have been studied: the homogeneous and the random history, respectively.The first one could be realized by temporarily subjecting the LC to a strong enough electric or magnetic field.The other case could be realized via a sudden quench from an isotropic ordering into the orientationally ordered phase.It has been demonstrated that for random histories  SRO is expected.On the other hand, for homogeneous histories LRO, QLRO, or even SRO could be realized depending on the values of  and .These observations are from qualitative and roughly also from quantitative point in line with known behavior of magnetic systems [24][25][26].

A. Calculation of Thermal Fluctuations
In order to set the thermal fluctuation vector Δ  ⇀   in (11), one must consider the rotation of the spins in their respective eigen frame denoted by unit vectors { ⃗  (0)  , ⃗  (0)  , ⃗  (0)  }.The axis is aligned along the spin direction.Only the rotations of a spin about the perpendicular local and -axes are relevant, and it is assumed that both rotations are independent of each other.We imagine the thermal fluctuations as the rotational kinetic free energy contributions with two independent components of angular velocity   and   .The related rotational kinetic energy reads where Ω  and Ω  are the corresponding rotation angles and  1 is the appropriate dimensionless tuning parameter.Assuming the canonical distribution we express the probability  for a range of angles within the phase-space element Ω  Ω  as where  is the normalization constant.It is convenient to rewrite (A.2) in the two-dimensional polar coordinate system.Taking into account the dimensionless quantities defined in Section 2 it follows where Ω is the magnitude of the vector ⃗ and the angle  defines its direction with respect to ⃗  (0)  .The distribution of  is uniform in the interval (0, 2), while the probability distribution function for Ω in the interval (0, ∝ ) is given by where  2 =  * / 1 .
After the thermal rotation of the spin in its own frame, the coordinates of the rotated spin are transformed into the laboratory system { ⃗   , ⃗   , ⃗   }.The thermal rotation of the individual spin can be finally written in the laboratory frame using the following vector transformation formula: The elements of the matrix are expressed with the spin components prior to rotations.Two rotations of the coordinate system are performed; the first for the spherical azimuthal angle  about the -axis, and the second one for the polar angle Ω about the new -axis.The vector to the right of the matrix is the rotated spin in its own coordinate system.

B. Impact of Impurities on I-N Phase Transition Temperature
We consider a mesomorphic binary mixture of rod-like molecules of comparable size [27,28].The concentration of the first component (i.e., LC spins) is given by 1 −  while the concentration of the second component (i.e., impurities) by .
The state of the th component is defined by the orientational distribution function   (Ω), where the angle Ω determines the orientation of a spin with respect to the symmetry axis.The low temperature orientationally ordered nematic phase is determined by an infinite set of orientational order parameters The dominant role is played by the order parameters  1 = (1) 2 and  2 = 2 , that is, where  11 is the intermolecular orientational interaction between components 1,  22 the orientational interaction between the components 2, while  12 is the orientational interaction between the spins of components 1 and 2, respectively.
The free energy difference between the nematic and isotropic phases can be written as where , , ,  * are positive material constants and   describes the contribution of impurities.In the temperature regime where the nematic phase is stable, the temperature dependence of constants , ,  can be neglected.Furthermore, if the interaction coupling with the impurities is negligible, then the I-N phase transition temperature is given by =  (0)  −  * , (B.7) where  (0)  determines the I-N phase transition for  = 0.

Figure 2 :
Figure 2: The hysteretic behavior of () is shown on increasing and decreasing  in different values of ,  = 1,  = 60.