# Hysteresis in Two-Dimensional Liquid Crystal Models

Authors: Slavko Buček; Samo Kralj; T. J. Sluckin
Journal: Advances in Condensed Matter Physics
 (2015)

---

## Abstract

We make a numerical study of hysteresis loop shapes within a generalized two-dimensional Random Anisotropy Nematic (RAN) model at zero temperature. The hysteresis loops appear on cycling a static external ordering field. Ordering in these systems is history dependent and involves interplay between the internal coupling constantJ, the anisotropy random field D, and the ordering external field H. Here the external field is represented by a traceless tensor, analogous to extension-type fields in continuum mechanics. The calculations use both a mean field approach and full lattice simulations. Our analysis suggests the existence of two qualitatively different solutions, which we denote as symmetric and symmetry breaking. For the set of parameters explored, only the symmetric solutions are stable. Both approaches yield qualitatively similar hysteresis curves, which are manifested either by single or double loops. But the quantitative differences indicate that mean field estimates are only of limited predictive value.


---

## Body

## 1. Introduction

Understanding the effect of disorder on structures reached via continuous symmetry breaking is of interest for various branches of physics [1, 2]. Such systems exhibit almost unavoidably topological defects [3, 4], which are in general stabilized by disorder [5–7]. The resulting structures are in general rich in metastable states, which are separated from each other by relatively high energy barriers with respect to thermal energies. Transitions between competing states often require local topological changes, which are in general energetically costly. Consequently, the sample history can be significant even on the macroscopic properties of the system.The pioneering investigations of spatially randomly perturbed system exhibiting continuous symmetry breaking were carried out in magnetism [8–10]. Other studies [11] show that disorder strongly influences phase behavior, as a result of the presence of Goldstone modes in gauge fields reflecting broken symmetries [11]. For example, the Larkin-Imry-Ma theorem [11, 12] claims that even infinitesimally weak random field-type disorder destabilizes long range order (LRO) with respect to short range order (SRO). However, later studies [13] demonstrate the occurrence rather of quasi-long-range-order (QLRO), characterized by algebraic decay of spatial correlations. A well-known example of such a system is the Bragg glass phase in dirty superconductors [14]. Glass-like properties are typically observed in these systems, which as a rule manifest history-dependent macroscopic responses when subject to an external ordering field [10]. Because of their potential applications in various memory storage devices, the hysteresis response [15, 16] to cycling external fields are of particular interest. Even in simple minimal models [17], the competition between the inherent ordering interactions, the disorder strength, and the external ordering field can lead to a relatively rich set of behaviors.In recent decades, there have been several studies of the influence of disorder on continuous symmetry broken phases in a number of different liquid crystal (LC) phases [6, 7, 18, 19]. These studies include a rich variety of different LC phases and also enable the type and strength of the disorder to be relatively well controlled. As such they represent a useful testing ground for the investigation of universal features related with disorder. We note that liquid crystals are relatively easily experimentally accessible due to their liquid character, softness, and optical anisotropy [20]. The simplest representative is the nematic uniaxial LC phase, whose local mesoscopic orientational ordering is typically represented by the nematic director field n, where the directions ±
n are physically equivalent [21, 22], and |
n
|
=
1. In the bulk, equilibrium nematic LCs tend to align parallel along a single symmetry direction, in analogous fashion to local magnetization m in ferromagnets. Random field-type disorder in LC structures can be imposed using a number of different procedures. Well-known examples are LCs confined to porous matrices such as controlled-pore glasses, aerogels or Russian glasses, or LC-aerosil mixtures [6, 18]. The latter systems are of particular interest because one can tune disorder strength, and even the type of disorder, by varying aerosil nanoparticle concentration [6, 23, 24].However, the symmetries ofn and m differ. Consequently, some features observed in magnetic and LC systems could also differ [19]. For example, magnetic systems exhibit line defects only with integer values of the corresponding winding number m [3, 20]. But the head-to-tail invariance of n permits nematic LCs also to exhibit line defects characterized by half-integer winding numbers.In this paper we consider the hysteresis properties of a randomly perturbed two-dimensional nematic LC described by the Random Anisotropy Nematic (RAN) [25] lattice model at zero temperature. The orientational order of the LC is measured by a tensor nematic order parameter. The external extensional ordering field is a traceless tensor and possesses the symmetry of an extensional stress of a type shown often in continuum mechanics [26]. We then calculate the global nematic response to the cycling of this field. Our particular interest lies in the existence of qualitatively different types of hysteresis loop. The calculations use both lattice-type simulations and also a mean field-type approach.The plan of the paper is as follows. In Section2 we formulate the problem. We introduce our two-dimensional lattice model in terms of tensor nematic order parameter at zero temperature. In Section 3 we derive mean field solution for average nematic ordering for a given external ordering field. Results are presented in Section 4 where we focus on possible shapes of hysteresis loops on cycling the static external field. We analyze both mean field behaviour and lattice model simulation results. In the final section we summarize the results of the calculations.


## 2. Formulation of Problem

### 2.1. Model

We consider ad
-dimensional lattice model of n dimensional nematic liquid crystal molecules. These are subject to external fields as well as local disordering random fields. The model is a generalization of the random anisotropy nematic (RAN) model discussed elsewhere by some of the present writers [25], which in turn is a random field generalization of the classic Lebwohl-Lasher lattice model [27] of a nematic liquid crystal. A rotor, whose direction is defined by n
≡
-
n, is placed at each site i. In this paper, d
=
n
=
2.Orientational properties of the rotors at sitei are determined by the traceless tensor quantity Q
i:
(1)
Q
α
β
i
=
n
α
i
n
β
i
-
1
n
δ
α
β
.
Greek indices label components of the Cartesian coordinate system, while Latin indices determine lattice sites. We will refer to the quantities Q
α
β
i astensor spins, to distinguish them from vector spins, whose hysteresis properties have been well studied in the literature [17].The Hamiltonian of the system reads(2)
H
=
-
1
2
J
∑
n
n
,
i
j
Q
α
β
i
Q
α
β
j
-
H
α
β
∑
i
Q
α
β
i
-
D
∑
i
N
α
β
i
Q
α
β
i
,
where repeated Greek indices in a product indicate the application of a summation convention.The meaning of the terms in (2) is as follows.(a)
The first term involves interactions between nearest neighbor sites and favours parallel orientation of neighboring LC molecules.J
>
0 is a measure of the interaction strength between the nearest neighbors (n
n). The factor 1
/
2 is introduced to avoid double counting. Later we will suppose z nearest neighbors, with z
=
4 in the square lattice, and z
=
6 in the hexagonal lattice. We refer to this term as theinternal field term.
(b)
The second term describes the coupling with an external tensor ordering fieldH
α
β. We define this in terms of a single parameter H, and a set of principal axes, which in this case coincide with x and y axes:
(3)
H
α
β
=
2
H
e
x
,
α
e
x
,
β
-
1
2
δ
α
β
=
-
2
H
e
y
,
α
e
y
,
β
-
1
2
δ
α
β
,
wheree
x is the unit vector in the x direction and e
y is the unit vector in the y direction. This tensor field takes the form of an extension field, of the type that frequently occurs in continuum mechanics. We address elsewhere in the paper the difference between this tensor form and the uniaxial random tensor field introduced by Cleaver et al. [25], but note at this stage only that for n
=
2 the two forms are completely equivalent. We will refer to this term as theexternal field term.
(c)
The final term represents the coupling between the local tensor spin and a time-independent quenched local field of uniform strengthD
>
0, where the tensor
(4)
N
α
β
i
=
e
α
i
e
β
i
-
1
2
δ
α
β
is expressed in terms of unit vectorse
i
≡
-
e
i. The directions e
i are drawn randomly from an isotropic orientational distribution, and the directions at different sites are uncorrelated. As in the case of the external field term, this tensor is an extensional field. We refer to this term as therandom field term. Note that the only difference between the terms in (3) and (4) is that the field in (3) aligns the spins in the same direction at each site, whereas in (4) the direction of the aligning field varies from site to site.


### 2.2. Angular Representation

In the(
x
,
y
) Cartesian frame, the dynamical state of the system can be parameterized at each site in terms of a local spin direction n
i
=
(
cos
⁡
⁡
θ
i
,
sin
⁡
θ
i
) and a locally favored orientation e
i
=
(
cos
⁡
⁡
α
i
,
sin
⁡
α
i
). Then it is readily shown that
(5)
Q
i
=
1
2
cos
⁡
⁡
2
θ
i
sin
⁡
2
θ
i
sin
⁡
2
θ
i
-
cos
⁡
⁡
2
θ
i
,
H
=
H
1
0
0
-
1
,
N
i
=
1
2
cos
⁡
⁡
2
α
i
sin
⁡
2
α
i
sin
⁡
2
α
i
cos
⁡
⁡
2
α
i
.
The energy can now be expressed in terms of the quantities θ
i
,
α
i. Substituting (5) into (2) we obtain
(6)
H
=
-
J
4
∑
n
n
,
i
j
cos
⁡
⁡
2
θ
i
-
θ
j
-
H
∑
i
cos
⁡
⁡
2
θ
i
-
D
2
∑
i
cos
⁡
⁡
2
θ
i
-
α
i
.
Although the model is originally expressed in terms of local tensor quantities, we note the similarity of this form to the classical random field X
Y vector spin model Hamiltonian
(7)
H
*
=
-
J
~
2
∑
n
n
,
i
j
cos
⁡
⁡
θ
i
-
θ
j
-
H
~
∑
i
cos
⁡
⁡
θ
i
-
D
~
∑
i
cos
⁡
⁡
(
θ
i
-
α
i
)
,
where now J
~ is the local spin interaction energy, H
~ is a magnetic field, and D
~ is the coupling with a random field direction. Apart from the factors of 2 in the angles and differences in convention concerning the values of parameters, (6) and (7) are extremely similar. We expect the results of statistical mechanism calculations to mirror this.We further observe that the field energy at a site is-
H
α
β
Q
α
β
i
=
-
H
cos
⁡
⁡
2
θ
i. Thus positive H prefers alignment in the x direction, while negative H prefers alignment in the y direction. Likewise the random energy at site i is given by -
D
N
α
β
i
Q
α
β
i
=
-
1
/
2
D
cos
⁡
⁡
2
θ
i
-
α
i, implying that for sufficiently large D,  θ
i
=
α
i.


### 2.3. Order Parameters

Our study concerns the average behaviour of the whole lattice. We can define a mean tensor order parameter(8)
S
=
Q
i
=
S
x
x
S
x
y
S
y
x
S
y
y
=
η
ζ
ζ
-
η
,
where ⋯ indicates spatial averaging. Thus the key order parameters are
(9a)
η
=
cos
⁡
⁡
2
θ
i
=
1
N
∑
i
cos
⁡
⁡
2
θ
i
,
(9b)
ζ
=
sin
⁡
2
θ
i
=
1
N
∑
i
sin
⁡
2
θ
i
,
where the averaging is over the all N sites. In the limits H
→
±
∞, the external field dominates the internal coupling and the random field terms, and hence η
H
→
+
∞
→
1, whereas η
H
→
-
∞
→
-
1. This corresponds to a rotation from the x to the y directions as H goes from +
∞ to -
∞.The quantityS responds to H, and thus we expect it to take the same form. The order parameter  η tends to adjust to the imposed field H. In this case the principal axes of S lie in the same directions as those of the forcing field H
α
β. By contrast, the quantity  ζ is an off-diagonal term in the principal frame of reference. In general, the symmetry of the system might suggest ζ
=
0, if we may suppose y and -
y to be equivalent. If ζ
≠
0, the symmetry of the system is broken, and we would in general expect a pair of solutions with positive and negative values. In this case, the principal axes of S will be rotated at an (plus or minus) oblique with respect to H. We henceforth refer to solutions with ζ
=
0 and ζ
≠
0 assymmetric solutions andsymmetry breaking solutions, respectively.


### 2.4. Procedure

The dynamics of the system atT
=
0 is governed by the instantaneous local field at site i, given by
(10)
G
α
β
i
=
-
δ
H
δ
Q
α
β
i
=
H
α
β
+
D
N
α
β
i
+
J
∑
n
n
,
j
≠
i
Q
α
β
j
,
where the local fields G
α
β
i are necessarily symmetric and traceless. Each individual tensor spin seeks to minimize the local energy
(11)
H
i
=
-
G
α
β
i
Q
α
β
i
,
where we bear in mind that ∑
i
H
i
≠
H. The sum of the local energies does not equal the total Hamiltonian, because of the need to include a compensating factor of +
1
/
2
J
∑
n
n
Q
α
β
i
Q
α
β
j, to avoid double counting of nearest neighbor interaction terms.The dynamics of the system proceeds in discrete time, such thatQ
α
β
i
t
+
1 minimizes the quantity -
G
α
β
i
t
Q
α
β
i
t
+
1. In the n
=
2, d
=
2 case discussed here, this only involves ensuring that(a)
the principal axes ofQ
α
β
i
t
+
1 coincide with those G
α
β
i
t;
(b)
the normalization ofQ
α
β
i
t
+
1, given by (1), is correct.We thus reach the result that
(12)
Q
i
t
+
1
=
1
2
cos
⁡
⁡
2
θ
i
t
+
1
sin
⁡
2
θ
i
t
+
1
sin
⁡
2
θ
i
t
+
1
-
cos
⁡
⁡
2
θ
i
t
+
1
=
1
tr
⁡
G
i
t
2
G
α
β
i
t
=
1
2
G
11
i
t
2
+
G
12
i
t
2
G
11
i
t
G
12
i
t
G
12
i
t
-
G
11
i
t
.
Given a lattice configuration of random tensor fields N
α
β
i and a starting configuration of spins Q
i
0, the discrete dynamics is iterated until a stable configuration, or equivalently a fixed point of the dynamics is reached. The condition for stability at site i, from (12), is then that
(13)
1
2
cos
⁡
⁡
2
θ
i
sin
⁡
2
θ
i
sin
⁡
2
θ
i
-
cos
⁡
⁡
2
θ
i
=
1
2
G
11
i
2
+
G
12
i
2
G
11
i
G
12
i
G
12
i
-
G
11
i
.
This corresponds to a local minimum of the energy, but not necessarily, of course, a global minimum. The local field combines the random field, the internal field, and the external field. Which of these dominates depends on a number of factors. If the magnitude of the external field is sufficiently large, other terms must be a small perturbation. Likewise, if the magnitude of the internal field is small compared to the random field, it can never play more than a minor role, even at low external fields. Finally, the history of the sample also plays a major role. Under some circumstances it is hard to destroy previously well-established order.We study the hysteretic behavior of ensembles of spins at zero temperature. In each numerical experiment we keep the same lattice with fixed value ofJ
>
0 and D
>
0 and fixed values of α
i at each site. The ordering field strength H is cycled between very large positive (H
=
H
max
⁡
>
0) and negative (H
=
-
H
max
⁡) values. The energy is first minimised at large positive H, corresponding to the vector e
i aligning along the x axis at all sites. The field is progressively reduced in small steps (the “down” part of the cycle). At each step, a new energy minimum is sought, starting from the previous equilibrium configuration, and following an energy-downhill procedure. Eventually, after many steps, H is large and negative, and all e
i essentially lie along the y direction, that is, rotated by π
/
2 with respect to their original direction. The process is then reversed, and H is taken from -
∞ to +
∞ (the “up” part of the cycle).We also develop mean field equations to describe this process and compare these results with the results of our simulations. Hysteresis corresponds to unambiguously different values of the order parameters during the up and down parts of the cycle. We also examine in detail the process of relaxation and tensor spin rotation.


## 2.1. Model

We consider ad
-dimensional lattice model of n dimensional nematic liquid crystal molecules. These are subject to external fields as well as local disordering random fields. The model is a generalization of the random anisotropy nematic (RAN) model discussed elsewhere by some of the present writers [25], which in turn is a random field generalization of the classic Lebwohl-Lasher lattice model [27] of a nematic liquid crystal. A rotor, whose direction is defined by n
≡
-
n, is placed at each site i. In this paper, d
=
n
=
2.Orientational properties of the rotors at sitei are determined by the traceless tensor quantity Q
i:
(1)
Q
α
β
i
=
n
α
i
n
β
i
-
1
n
δ
α
β
.
Greek indices label components of the Cartesian coordinate system, while Latin indices determine lattice sites. We will refer to the quantities Q
α
β
i astensor spins, to distinguish them from vector spins, whose hysteresis properties have been well studied in the literature [17].The Hamiltonian of the system reads(2)
H
=
-
1
2
J
∑
n
n
,
i
j
Q
α
β
i
Q
α
β
j
-
H
α
β
∑
i
Q
α
β
i
-
D
∑
i
N
α
β
i
Q
α
β
i
,
where repeated Greek indices in a product indicate the application of a summation convention.The meaning of the terms in (2) is as follows.(a)
The first term involves interactions between nearest neighbor sites and favours parallel orientation of neighboring LC molecules.J
>
0 is a measure of the interaction strength between the nearest neighbors (n
n). The factor 1
/
2 is introduced to avoid double counting. Later we will suppose z nearest neighbors, with z
=
4 in the square lattice, and z
=
6 in the hexagonal lattice. We refer to this term as theinternal field term.
(b)
The second term describes the coupling with an external tensor ordering fieldH
α
β. We define this in terms of a single parameter H, and a set of principal axes, which in this case coincide with x and y axes:
(3)
H
α
β
=
2
H
e
x
,
α
e
x
,
β
-
1
2
δ
α
β
=
-
2
H
e
y
,
α
e
y
,
β
-
1
2
δ
α
β
,
wheree
x is the unit vector in the x direction and e
y is the unit vector in the y direction. This tensor field takes the form of an extension field, of the type that frequently occurs in continuum mechanics. We address elsewhere in the paper the difference between this tensor form and the uniaxial random tensor field introduced by Cleaver et al. [25], but note at this stage only that for n
=
2 the two forms are completely equivalent. We will refer to this term as theexternal field term.
(c)
The final term represents the coupling between the local tensor spin and a time-independent quenched local field of uniform strengthD
>
0, where the tensor
(4)
N
α
β
i
=
e
α
i
e
β
i
-
1
2
δ
α
β
is expressed in terms of unit vectorse
i
≡
-
e
i. The directions e
i are drawn randomly from an isotropic orientational distribution, and the directions at different sites are uncorrelated. As in the case of the external field term, this tensor is an extensional field. We refer to this term as therandom field term. Note that the only difference between the terms in (3) and (4) is that the field in (3) aligns the spins in the same direction at each site, whereas in (4) the direction of the aligning field varies from site to site.


## 2.2. Angular Representation

In the(
x
,
y
) Cartesian frame, the dynamical state of the system can be parameterized at each site in terms of a local spin direction n
i
=
(
cos
⁡
⁡
θ
i
,
sin
⁡
θ
i
) and a locally favored orientation e
i
=
(
cos
⁡
⁡
α
i
,
sin
⁡
α
i
). Then it is readily shown that
(5)
Q
i
=
1
2
cos
⁡
⁡
2
θ
i
sin
⁡
2
θ
i
sin
⁡
2
θ
i
-
cos
⁡
⁡
2
θ
i
,
H
=
H
1
0
0
-
1
,
N
i
=
1
2
cos
⁡
⁡
2
α
i
sin
⁡
2
α
i
sin
⁡
2
α
i
cos
⁡
⁡
2
α
i
.
The energy can now be expressed in terms of the quantities θ
i
,
α
i. Substituting (5) into (2) we obtain
(6)
H
=
-
J
4
∑
n
n
,
i
j
cos
⁡
⁡
2
θ
i
-
θ
j
-
H
∑
i
cos
⁡
⁡
2
θ
i
-
D
2
∑
i
cos
⁡
⁡
2
θ
i
-
α
i
.
Although the model is originally expressed in terms of local tensor quantities, we note the similarity of this form to the classical random field X
Y vector spin model Hamiltonian
(7)
H
*
=
-
J
~
2
∑
n
n
,
i
j
cos
⁡
⁡
θ
i
-
θ
j
-
H
~
∑
i
cos
⁡
⁡
θ
i
-
D
~
∑
i
cos
⁡
⁡
(
θ
i
-
α
i
)
,
where now J
~ is the local spin interaction energy, H
~ is a magnetic field, and D
~ is the coupling with a random field direction. Apart from the factors of 2 in the angles and differences in convention concerning the values of parameters, (6) and (7) are extremely similar. We expect the results of statistical mechanism calculations to mirror this.We further observe that the field energy at a site is-
H
α
β
Q
α
β
i
=
-
H
cos
⁡
⁡
2
θ
i. Thus positive H prefers alignment in the x direction, while negative H prefers alignment in the y direction. Likewise the random energy at site i is given by -
D
N
α
β
i
Q
α
β
i
=
-
1
/
2
D
cos
⁡
⁡
2
θ
i
-
α
i, implying that for sufficiently large D,  θ
i
=
α
i.


## 2.3. Order Parameters

Our study concerns the average behaviour of the whole lattice. We can define a mean tensor order parameter(8)
S
=
Q
i
=
S
x
x
S
x
y
S
y
x
S
y
y
=
η
ζ
ζ
-
η
,
where ⋯ indicates spatial averaging. Thus the key order parameters are
(9a)
η
=
cos
⁡
⁡
2
θ
i
=
1
N
∑
i
cos
⁡
⁡
2
θ
i
,
(9b)
ζ
=
sin
⁡
2
θ
i
=
1
N
∑
i
sin
⁡
2
θ
i
,
where the averaging is over the all N sites. In the limits H
→
±
∞, the external field dominates the internal coupling and the random field terms, and hence η
H
→
+
∞
→
1, whereas η
H
→
-
∞
→
-
1. This corresponds to a rotation from the x to the y directions as H goes from +
∞ to -
∞.The quantityS responds to H, and thus we expect it to take the same form. The order parameter  η tends to adjust to the imposed field H. In this case the principal axes of S lie in the same directions as those of the forcing field H
α
β. By contrast, the quantity  ζ is an off-diagonal term in the principal frame of reference. In general, the symmetry of the system might suggest ζ
=
0, if we may suppose y and -
y to be equivalent. If ζ
≠
0, the symmetry of the system is broken, and we would in general expect a pair of solutions with positive and negative values. In this case, the principal axes of S will be rotated at an (plus or minus) oblique with respect to H. We henceforth refer to solutions with ζ
=
0 and ζ
≠
0 assymmetric solutions andsymmetry breaking solutions, respectively.


## 2.4. Procedure

The dynamics of the system atT
=
0 is governed by the instantaneous local field at site i, given by
(10)
G
α
β
i
=
-
δ
H
δ
Q
α
β
i
=
H
α
β
+
D
N
α
β
i
+
J
∑
n
n
,
j
≠
i
Q
α
β
j
,
where the local fields G
α
β
i are necessarily symmetric and traceless. Each individual tensor spin seeks to minimize the local energy
(11)
H
i
=
-
G
α
β
i
Q
α
β
i
,
where we bear in mind that ∑
i
H
i
≠
H. The sum of the local energies does not equal the total Hamiltonian, because of the need to include a compensating factor of +
1
/
2
J
∑
n
n
Q
α
β
i
Q
α
β
j, to avoid double counting of nearest neighbor interaction terms.The dynamics of the system proceeds in discrete time, such thatQ
α
β
i
t
+
1 minimizes the quantity -
G
α
β
i
t
Q
α
β
i
t
+
1. In the n
=
2, d
=
2 case discussed here, this only involves ensuring that(a)
the principal axes ofQ
α
β
i
t
+
1 coincide with those G
α
β
i
t;
(b)
the normalization ofQ
α
β
i
t
+
1, given by (1), is correct.We thus reach the result that
(12)
Q
i
t
+
1
=
1
2
cos
⁡
⁡
2
θ
i
t
+
1
sin
⁡
2
θ
i
t
+
1
sin
⁡
2
θ
i
t
+
1
-
cos
⁡
⁡
2
θ
i
t
+
1
=
1
tr
⁡
G
i
t
2
G
α
β
i
t
=
1
2
G
11
i
t
2
+
G
12
i
t
2
G
11
i
t
G
12
i
t
G
12
i
t
-
G
11
i
t
.
Given a lattice configuration of random tensor fields N
α
β
i and a starting configuration of spins Q
i
0, the discrete dynamics is iterated until a stable configuration, or equivalently a fixed point of the dynamics is reached. The condition for stability at site i, from (12), is then that
(13)
1
2
cos
⁡
⁡
2
θ
i
sin
⁡
2
θ
i
sin
⁡
2
θ
i
-
cos
⁡
⁡
2
θ
i
=
1
2
G
11
i
2
+
G
12
i
2
G
11
i
G
12
i
G
12
i
-
G
11
i
.
This corresponds to a local minimum of the energy, but not necessarily, of course, a global minimum. The local field combines the random field, the internal field, and the external field. Which of these dominates depends on a number of factors. If the magnitude of the external field is sufficiently large, other terms must be a small perturbation. Likewise, if the magnitude of the internal field is small compared to the random field, it can never play more than a minor role, even at low external fields. Finally, the history of the sample also plays a major role. Under some circumstances it is hard to destroy previously well-established order.We study the hysteretic behavior of ensembles of spins at zero temperature. In each numerical experiment we keep the same lattice with fixed value ofJ
>
0 and D
>
0 and fixed values of α
i at each site. The ordering field strength H is cycled between very large positive (H
=
H
max
⁡
>
0) and negative (H
=
-
H
max
⁡) values. The energy is first minimised at large positive H, corresponding to the vector e
i aligning along the x axis at all sites. The field is progressively reduced in small steps (the “down” part of the cycle). At each step, a new energy minimum is sought, starting from the previous equilibrium configuration, and following an energy-downhill procedure. Eventually, after many steps, H is large and negative, and all e
i essentially lie along the y direction, that is, rotated by π
/
2 with respect to their original direction. The process is then reversed, and H is taken from -
∞ to +
∞ (the “up” part of the cycle).We also develop mean field equations to describe this process and compare these results with the results of our simulations. Hysteresis corresponds to unambiguously different values of the order parameters during the up and down parts of the cycle. We also examine in detail the process of relaxation and tensor spin rotation.


## 3. Mean Field Approximation

### 3.1. The Local Field

In the following we calculate ordering properties of an ensemble driven by Hamiltonian defined in (2), using a mean field approximation valid in the limit that all sites interact equally with all others. This approach follows previous work by Strogatz et al. [28], who were interested in conduction in charge density waves, and, by Shukla and Kharwanlang [17], who focused on the much more closely related problem of random vector spins.A mean field theory requires the determination of self-consistent local field. This involves ignoring correlations between sites in (10), leading to a mean local field at site i given by
(14)
G
-
α
β
i
=
H
α
β
+
D
N
α
β
i
+
J
z
S
α
β
.


### 3.2. Symmetric Solution

Assuming no symmetries are broken and taking into accountS
x
x
+
S
y
y
=
0 we can define S using a single variable. From the definition (9a) it follows that
(15)
S
=
1
2
η
0
0
-
η
,
G
=
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
1
2
D
sin
⁡
2
α
i
1
2
D
sin
⁡
2
α
i
-
H
-
1
2
D
cos
⁡
⁡
2
α
i
-
1
2
J
z
η
.
Within a mean field philosophy, we may suppose that Q
α
β
i is parallel to G
α
β
i at each site. Hence
(16)
Q
α
β
i
=
1
2
cos
⁡
⁡
2
θ
i
sin
⁡
2
θ
i
sin
⁡
2
θ
i
-
cos
⁡
⁡
2
θ
i
=
λ
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
1
2
D
sin
⁡
2
α
i
1
2
D
sin
⁡
2
α
i
-
H
-
1
2
D
cos
⁡
⁡
2
α
i
-
1
2
J
z
η
,
and hence the equation for θ
i is
(17)
cos
⁡
2
θ
i
=
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
·
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
2
+
1
4
D
2
sin
2
⁡
2
α
i
-
1
/
2
.
The reader should note the strong analogy between procedure outlined above and other work in the mean field theory of disordered systems at zero temperature (see, e.g., [17, 28]) The philosophy is similar to that of mean field theory in temperature-dependent nonglassy materials. However, unlike in the case of temperature-dependent nonrandom systems, as far as we are aware, there is no detailed theory discussing the accuracy and dimensional dependence of mean field theories in these disordered systems. There is, of course, an expectation, based on the applicability of mean field theory in temperature-dependent nonglassy materials, that mean field theory becomes more accurate in the limits of long-ranged interactions and high dimensionality. At this stage, we merely present the mean field results because the calculation is heuristically appealing and at this stage the only approximation available.The self-consistent equation forη involves taking the average, which we do by integration:
(18)
η
=
1
π
∫
0
π
d
α
cos
⁡
⁡
2
θ
i
α
,
leading to a final self-consistent equation of
(19)
η
=
1
π
∫
0
π
d
α
·
H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
2
+
1
4
D
2
sin
2
⁡
2
α
-
1
/
2

H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
lllllllllllllllll
·
H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
2
llllllllllllllllllll
+
1
4
D
2
sin
2
⁡
2
α
-
1
/
2
.
This equation is exactly of the same form as that in the Shukla-Kharwanlang X
Y model solution [17], apart from scaling.


### 3.3. Symmetry Breaking Solution

Note that we might expect the liquid crystal direction to swing around, rather than just to change fromx to y, so that the whole system swings around coherently. Then we need to take into account also the order parameter  ζ; see (9b). Now (8) becomes
(20)
S
=
η
ζ
ζ
-
η
,
G
=
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
1
2
D
sin
⁡
2
α
i
+
1
2
J
z
ζ
1
2
D
sin
⁡
2
α
i
+
1
2
J
z
ζ
-
H
-
1
2
D
cos
⁡
⁡
2
α
i
-
1
2
J
z
η
.
Requesting parallel orientation of Q
i and G it follows
(21)
Q
α
β
i
=
1
2
cos
⁡
⁡
2
θ
i
sin
⁡
2
θ
i
sin
⁡
2
θ
i
-
cos
⁡
⁡
2
θ
i
=
λ
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
1
2
D
sin
⁡
2
α
i
+
1
2
J
z
ζ
1
2
D
sin
⁡
2
α
i
+
1
2
J
z
ζ
-
H
-
1
2
D
cos
⁡
⁡
2
α
i
-
1
2
J
z
η
.
The (two) self-consistent equation is now
(22)
η
=
1
π
∫
0
π
d
α
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
2
-
1
/
2

H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
lllllllllllllllll
·
H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
2
llllllllllllllllllll
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
2
-
1
/
2
,
ζ
=
1
π
∫
0
π
d
α
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
2
-
1
/
2

H
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
lllllllllllllllll
·
H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
2
llllllllllllllllllll
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
2
-
1
/
2
.


### 3.4. Effective Energy

The self-consistent equations can be derived from an effective energy per site which we define as(23)
ε
η
,
ζ
=
1
4
J
z
η
2
+
ζ
2
-
E
dis
η
,
ζ
,
E
dis
=
1
π
·
∫
0
π
d
α
H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
2
llllllllllllll
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
2
1
/
2
.
With this definition of the effective energy, the self-consistent equations become
(24)
∂
ε
∂
η
=
0
;
∂
ε
∂
ζ
=
0
.
In the restricted theory (ζ
=
0), corresponding tosymmetric solutions, the coercive field occurs when
(25)
∂
2
ε
∂
η
2
=
0
.


## 3.1. The Local Field

In the following we calculate ordering properties of an ensemble driven by Hamiltonian defined in (2), using a mean field approximation valid in the limit that all sites interact equally with all others. This approach follows previous work by Strogatz et al. [28], who were interested in conduction in charge density waves, and, by Shukla and Kharwanlang [17], who focused on the much more closely related problem of random vector spins.A mean field theory requires the determination of self-consistent local field. This involves ignoring correlations between sites in (10), leading to a mean local field at site i given by
(14)
G
-
α
β
i
=
H
α
β
+
D
N
α
β
i
+
J
z
S
α
β
.


## 3.2. Symmetric Solution

Assuming no symmetries are broken and taking into accountS
x
x
+
S
y
y
=
0 we can define S using a single variable. From the definition (9a) it follows that
(15)
S
=
1
2
η
0
0
-
η
,
G
=
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
1
2
D
sin
⁡
2
α
i
1
2
D
sin
⁡
2
α
i
-
H
-
1
2
D
cos
⁡
⁡
2
α
i
-
1
2
J
z
η
.
Within a mean field philosophy, we may suppose that Q
α
β
i is parallel to G
α
β
i at each site. Hence
(16)
Q
α
β
i
=
1
2
cos
⁡
⁡
2
θ
i
sin
⁡
2
θ
i
sin
⁡
2
θ
i
-
cos
⁡
⁡
2
θ
i
=
λ
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
1
2
D
sin
⁡
2
α
i
1
2
D
sin
⁡
2
α
i
-
H
-
1
2
D
cos
⁡
⁡
2
α
i
-
1
2
J
z
η
,
and hence the equation for θ
i is
(17)
cos
⁡
2
θ
i
=
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
·
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
2
+
1
4
D
2
sin
2
⁡
2
α
i
-
1
/
2
.
The reader should note the strong analogy between procedure outlined above and other work in the mean field theory of disordered systems at zero temperature (see, e.g., [17, 28]) The philosophy is similar to that of mean field theory in temperature-dependent nonglassy materials. However, unlike in the case of temperature-dependent nonrandom systems, as far as we are aware, there is no detailed theory discussing the accuracy and dimensional dependence of mean field theories in these disordered systems. There is, of course, an expectation, based on the applicability of mean field theory in temperature-dependent nonglassy materials, that mean field theory becomes more accurate in the limits of long-ranged interactions and high dimensionality. At this stage, we merely present the mean field results because the calculation is heuristically appealing and at this stage the only approximation available.The self-consistent equation forη involves taking the average, which we do by integration:
(18)
η
=
1
π
∫
0
π
d
α
cos
⁡
⁡
2
θ
i
α
,
leading to a final self-consistent equation of
(19)
η
=
1
π
∫
0
π
d
α
·
H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
2
+
1
4
D
2
sin
2
⁡
2
α
-
1
/
2

H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
lllllllllllllllll
·
H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
2
llllllllllllllllllll
+
1
4
D
2
sin
2
⁡
2
α
-
1
/
2
.
This equation is exactly of the same form as that in the Shukla-Kharwanlang X
Y model solution [17], apart from scaling.


## 3.3. Symmetry Breaking Solution

Note that we might expect the liquid crystal direction to swing around, rather than just to change fromx to y, so that the whole system swings around coherently. Then we need to take into account also the order parameter  ζ; see (9b). Now (8) becomes
(20)
S
=
η
ζ
ζ
-
η
,
G
=
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
1
2
D
sin
⁡
2
α
i
+
1
2
J
z
ζ
1
2
D
sin
⁡
2
α
i
+
1
2
J
z
ζ
-
H
-
1
2
D
cos
⁡
⁡
2
α
i
-
1
2
J
z
η
.
Requesting parallel orientation of Q
i and G it follows
(21)
Q
α
β
i
=
1
2
cos
⁡
⁡
2
θ
i
sin
⁡
2
θ
i
sin
⁡
2
θ
i
-
cos
⁡
⁡
2
θ
i
=
λ
H
+
1
2
D
cos
⁡
⁡
2
α
i
+
1
2
J
z
η
1
2
D
sin
⁡
2
α
i
+
1
2
J
z
ζ
1
2
D
sin
⁡
2
α
i
+
1
2
J
z
ζ
-
H
-
1
2
D
cos
⁡
⁡
2
α
i
-
1
2
J
z
η
.
The (two) self-consistent equation is now
(22)
η
=
1
π
∫
0
π
d
α
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
2
-
1
/
2

H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
lllllllllllllllll
·
H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
2
llllllllllllllllllll
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
2
-
1
/
2
,
ζ
=
1
π
∫
0
π
d
α
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
2
-
1
/
2

H
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
lllllllllllllllll
·
H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
2
llllllllllllllllllll
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
2
-
1
/
2
.


## 3.4. Effective Energy

The self-consistent equations can be derived from an effective energy per site which we define as(23)
ε
η
,
ζ
=
1
4
J
z
η
2
+
ζ
2
-
E
dis
η
,
ζ
,
E
dis
=
1
π
·
∫
0
π
d
α
H
+
1
2
D
cos
⁡
⁡
2
α
+
1
2
J
z
η
2
llllllllllllll
+
1
2
D
sin
⁡
2
α
+
1
2
J
z
ζ
2
1
/
2
.
With this definition of the effective energy, the self-consistent equations become
(24)
∂
ε
∂
η
=
0
;
∂
ε
∂
ζ
=
0
.
In the restricted theory (ζ
=
0), corresponding tosymmetric solutions, the coercive field occurs when
(25)
∂
2
ε
∂
η
2
=
0
.


## 4. Numerical Results

We fix the value ofJ
>
0 and D
>
0 at zero temperature and focus on the hysteretic behavior of ensembles on cycling the ordering field strength H. We first analyze behavior using the mean field approximation. Afterwards we perform a more accurate lattice model simulation.We determine the ordering properties of the system on cyclingH from +
H
max
⁡
>
0 to -
H
max
⁡ and back to +
H
max
⁡. We solve equations iteratively for a given set of external parameters. A simulation is stopped when the difference between the new and old configuration is sufficiently small. In this context the criterion of smallness is that at each site the absolute difference between the previous and the updated configuration is below 10−5.In practice we begin simulations at a high enough value ofH for which η
~
1 and ζ
=
0. Next we decrease the field for a small enough step Δ
H and calculate the new solution by using the previous solution as an initial guess. In such a way we obtain {
η
(
H
)
,
ζ
(
H
)
} curves by carrying out one cycle in H values. We choose a large enough step Δ
H which does not affect {
η
(
H
)
,
ζ
(
H
)
} profiles.


### 4.1. Mean Field Results

We solve (22) iteratively using the Mathematica software [29]. For a specific set of parameters we stop iterations when the differences η
i
+
1
-
η
i and ζ
i
+
1
-
ζ
i between the ith and i
+
1th iteration steps are below 10−5.Qualitatively different solutions on cyclingH are shown in Figure 1. In Figure 1(a) we vary D for J
=
1 and conversely in Figure 1(b) we vary J for D
=
1. In all cases studied we obtain onlysymmetric solutions; that is, we find that ζ
=
0 within numerical accuracy for all set of parameters and system histories. In dependence of the ratio μ
=
D
/
J we observe either (i) no hysteresis, (ii) a single hysteresis loop, or (iii) a double hysteresis loop, in η
(
H
). For relatively low values  μ we obtain hysteretic behavior, where values of  η switch at a critical value of H
=
H
c
(
μ
) between η
~
±
1. With increasing values of  μ departures from a rectangularly shaped hysteretic profile become increasingly more pronounced. In addition the width of the hysteresis loop decreases. On further increasing  μ, a double hysteretic loop begins to occur at μ
=
μ
c
(
1
) (e.g., μ
c
(
1
)
=
D
c
1
/
J
~
0.7 for J
=
1). This persists till μ
=
μ
c
(
2
) (e.g., μ
c
(
2
)
=
D
c
2
/
J
~
1.2 for J
=
1), beyond which we observe gradual evolution in η
(
H
) without hysteresis.Mean field approximation. Dependence ofη
(
H
) on external field H through a complete hysteresis cycle. Either no hysteresis, single loop hysteresis, or double loop hysteresis is observed. (a) J
=
1; dashed line: D
=
0.3, full line: D
=
0.8, and dotted line: D
=
1.2. (b) D
=
1; full line: J
=
0.5, dashed line: J
=
1.2, and dotted line: J
=
2.
(a)
(b)In Figure2 we plot the effective energy per site  ε versus  η for specific values of D and H; J
=
1. The energy is calculated using (23). In Figure 2(a) we demonstrate the case corresponding to the coercive field H
c defined by (25). In the case shown H
c = 0.22 (H
c
=
-
0.22) for the “up” (down) part of the cycle for D
=
0.3. The corresponding hysteresis loop is plotted in Figure 1(a) with the dashed line. For H
<
H
c the ε
(
η
) dependence exhibits double minima and at H
=
H
c one minimum ceases to exist. In Figure 2(b) we show typical free energy landscape in the regime where double hysteresis loop exists. In the case shown D
=
0.8 and H
=
±
0.09 corresponding to the double hysteresis loop depicted with the full line in Figure 1(a).The mean field effective energy per siteε versus  η for different values of H and D, J
=
1,symmetric solutions. (a) D
=
0.3; full line: H
=
H
c
=
0.22 obtained for the “up” part of the single loop hysteresis cycle; dashed line: H
c
=
-
0.22 obtained for the “down” part of the cycle. Note that the physically sensible regime of  η is restricted to the interval [
-
1,1
]. (b) D
=
0.8; full line: H
=
0.09 obtained for the “up” part of the double loop hysteresis cycle; dashed line: H
=
-
0.09 obtained for the “down” part of the cycle.
(a)
(b)The observed behavior can be qualitatively understood in the following way. In the limitμ
≫
1 the disorder dominantly influences structural behavior if H
<
D. Therefore, for H
=
0 it is always the case that η
=
0 for a large enough system because the disorder enforces isotropic symmetry. A finite value of H breaks the symmetry of the system and consequently η
≠
0. Because J
≪
D some degree of local disorder field preference prevails in the regime D
<
H. Consequently the collective ordering tendency favored by the interaction constant J is always overwhelmed, either by the external field H or by the local random nematic disorder D. Conversely, in the limit μ
~
0 the collective behavior tendency favored by the constant J is large if J
<
H. Therefore, for a relatively weak ordering field H the system tends to be aligned along a symmetry breaking direction, which is history dependent, giving rise to a pronounced single loop hysteresis behavior.


### 4.2. Lattice Simulation

In the lattice simulation we calculate the local field tensor defined by (10) and align Q
α
β
i parallel to it. In one sweep we update in such a way all the lattice sites of a system. The sweeps are repeated until at each site the difference tr
⁡
Q
α
β
i
,
new
-
Q
α
β
i
,
old between the previous Q
α
β
i
,
old and new Q
α
β
i
,
new configuration is below 10−5, signalling that a fixed point structure is reached. From it we calculate order parameters  η and ζ using (9a) and (9b). We carried simulations for the system size 316
×
316 (≈105) for which the finite size effects were negligible.Simulations confirmed our MFA results regarding the absence of thesymmetry breaking solutions. Furthermore, also in this case three topologically different responses on cycling H were obtained (i.e., no hysteresis, single hysteresis loop, or double hysteresis loops). However, hysteresis behavior on varying the ratio μ
=
D
/
J exhibits quantitatively and in some regimes also qualitatively different behavior.Typical behavior on varyingμ
=
D
/
J is shown in Figure 3, where field cycles are calculated for D
=
1. On increasing J (i.e., decreasing  μ) we observe the following behavior. For  μ
≪
1 we obtain roughly rectangularly shaped hysteresis loop, where configurations switch between μ
=
±
1. Initially on decreasing  μ the hysteresis loop width Δ
H is gradually shrinking towards zero. Within numerical accuracy we obtain Δ
H
=
0 in the interval between μ
~
0.01 and μ
=
μ
c
(
2
)
~
0.1. For μ
>
μ
c
(
2
) a double hysteresis loop appears which persists till μ
c
(
1
)
~
0.22. Above this value we obtain S-shaped hysteresis loop where Δ
H gradually decreases. Finally, above μ
c
(
0
)
~
3.3 no hysteresis is observed.Figure 3
Lattice simulationD
=
1: dependence of η
(
H
) on external field H through a complete hysteresis cycle. We observe either (a) μ
=
4 (dashed line), no hysteresis; (b) μ
=
0.05 (rectangularly shaped), single loop hysteresis; (c) μ
=
0.15 (dash-dotted line), double loop hysteresis; or (d) μ
=
2, S-shaped single loop hysteresis.


## 4.1. Mean Field Results

We solve (22) iteratively using the Mathematica software [29]. For a specific set of parameters we stop iterations when the differences η
i
+
1
-
η
i and ζ
i
+
1
-
ζ
i between the ith and i
+
1th iteration steps are below 10−5.Qualitatively different solutions on cyclingH are shown in Figure 1. In Figure 1(a) we vary D for J
=
1 and conversely in Figure 1(b) we vary J for D
=
1. In all cases studied we obtain onlysymmetric solutions; that is, we find that ζ
=
0 within numerical accuracy for all set of parameters and system histories. In dependence of the ratio μ
=
D
/
J we observe either (i) no hysteresis, (ii) a single hysteresis loop, or (iii) a double hysteresis loop, in η
(
H
). For relatively low values  μ we obtain hysteretic behavior, where values of  η switch at a critical value of H
=
H
c
(
μ
) between η
~
±
1. With increasing values of  μ departures from a rectangularly shaped hysteretic profile become increasingly more pronounced. In addition the width of the hysteresis loop decreases. On further increasing  μ, a double hysteretic loop begins to occur at μ
=
μ
c
(
1
) (e.g., μ
c
(
1
)
=
D
c
1
/
J
~
0.7 for J
=
1). This persists till μ
=
μ
c
(
2
) (e.g., μ
c
(
2
)
=
D
c
2
/
J
~
1.2 for J
=
1), beyond which we observe gradual evolution in η
(
H
) without hysteresis.Mean field approximation. Dependence ofη
(
H
) on external field H through a complete hysteresis cycle. Either no hysteresis, single loop hysteresis, or double loop hysteresis is observed. (a) J
=
1; dashed line: D
=
0.3, full line: D
=
0.8, and dotted line: D
=
1.2. (b) D
=
1; full line: J
=
0.5, dashed line: J
=
1.2, and dotted line: J
=
2.
(a)
(b)In Figure2 we plot the effective energy per site  ε versus  η for specific values of D and H; J
=
1. The energy is calculated using (23). In Figure 2(a) we demonstrate the case corresponding to the coercive field H
c defined by (25). In the case shown H
c = 0.22 (H
c
=
-
0.22) for the “up” (down) part of the cycle for D
=
0.3. The corresponding hysteresis loop is plotted in Figure 1(a) with the dashed line. For H
<
H
c the ε
(
η
) dependence exhibits double minima and at H
=
H
c one minimum ceases to exist. In Figure 2(b) we show typical free energy landscape in the regime where double hysteresis loop exists. In the case shown D
=
0.8 and H
=
±
0.09 corresponding to the double hysteresis loop depicted with the full line in Figure 1(a).The mean field effective energy per siteε versus  η for different values of H and D, J
=
1,symmetric solutions. (a) D
=
0.3; full line: H
=
H
c
=
0.22 obtained for the “up” part of the single loop hysteresis cycle; dashed line: H
c
=
-
0.22 obtained for the “down” part of the cycle. Note that the physically sensible regime of  η is restricted to the interval [
-
1,1
]. (b) D
=
0.8; full line: H
=
0.09 obtained for the “up” part of the double loop hysteresis cycle; dashed line: H
=
-
0.09 obtained for the “down” part of the cycle.
(a)
(b)The observed behavior can be qualitatively understood in the following way. In the limitμ
≫
1 the disorder dominantly influences structural behavior if H
<
D. Therefore, for H
=
0 it is always the case that η
=
0 for a large enough system because the disorder enforces isotropic symmetry. A finite value of H breaks the symmetry of the system and consequently η
≠
0. Because J
≪
D some degree of local disorder field preference prevails in the regime D
<
H. Consequently the collective ordering tendency favored by the interaction constant J is always overwhelmed, either by the external field H or by the local random nematic disorder D. Conversely, in the limit μ
~
0 the collective behavior tendency favored by the constant J is large if J
<
H. Therefore, for a relatively weak ordering field H the system tends to be aligned along a symmetry breaking direction, which is history dependent, giving rise to a pronounced single loop hysteresis behavior.


## 4.2. Lattice Simulation

In the lattice simulation we calculate the local field tensor defined by (10) and align Q
α
β
i parallel to it. In one sweep we update in such a way all the lattice sites of a system. The sweeps are repeated until at each site the difference tr
⁡
Q
α
β
i
,
new
-
Q
α
β
i
,
old between the previous Q
α
β
i
,
old and new Q
α
β
i
,
new configuration is below 10−5, signalling that a fixed point structure is reached. From it we calculate order parameters  η and ζ using (9a) and (9b). We carried simulations for the system size 316
×
316 (≈105) for which the finite size effects were negligible.Simulations confirmed our MFA results regarding the absence of thesymmetry breaking solutions. Furthermore, also in this case three topologically different responses on cycling H were obtained (i.e., no hysteresis, single hysteresis loop, or double hysteresis loops). However, hysteresis behavior on varying the ratio μ
=
D
/
J exhibits quantitatively and in some regimes also qualitatively different behavior.Typical behavior on varyingμ
=
D
/
J is shown in Figure 3, where field cycles are calculated for D
=
1. On increasing J (i.e., decreasing  μ) we observe the following behavior. For  μ
≪
1 we obtain roughly rectangularly shaped hysteresis loop, where configurations switch between μ
=
±
1. Initially on decreasing  μ the hysteresis loop width Δ
H is gradually shrinking towards zero. Within numerical accuracy we obtain Δ
H
=
0 in the interval between μ
~
0.01 and μ
=
μ
c
(
2
)
~
0.1. For μ
>
μ
c
(
2
) a double hysteresis loop appears which persists till μ
c
(
1
)
~
0.22. Above this value we obtain S-shaped hysteresis loop where Δ
H gradually decreases. Finally, above μ
c
(
0
)
~
3.3 no hysteresis is observed.Figure 3
Lattice simulationD
=
1: dependence of η
(
H
) on external field H through a complete hysteresis cycle. We observe either (a) μ
=
4 (dashed line), no hysteresis; (b) μ
=
0.05 (rectangularly shaped), single loop hysteresis; (c) μ
=
0.15 (dash-dotted line), double loop hysteresis; or (d) μ
=
2, S-shaped single loop hysteresis.


## 5. Conclusion

In our numerical study we consider two-dimensional lattice model of nematic LC ordering in presence of an external ordering field and local disordering random field. We use a generalized Random Anisotropy Nematic model where we introduce an external tensor ordering field of the type frequently encountered in continuum mechanics. Our approach resembles work by Shukla and Kharwanlang [17], who examined the shape of hysteresis loops on cycling a statical external ordering field using random field X
Y model at zero temperature. Despite different symmetries of order parameters the latter approach and our model yield similar Hamiltonians in angular presentation of orientational ordering.We first calculated ordering properties of generalized RAN model using the mean field approach. We derived a self-consistent expression for average scalar nematic order parameters. In general we obtain two qualitatively different solutions to which we refer assymmetric solutionsandbroken symmetry solutions. In the first type of solutions the principal axes of mean nematic tensor order parameter and of the external tensor field coincide. On contrary in the latter type of solutions the system’s nematic order parameter possesses off diagonal terms within the principal frame of reference imposed by H. For the set of parameters explored in our simulations we obtained onlysymmetric solutions, well represented by single scalar order parameter  η.We were particularly interested in hysteresis curves resulting when the forcing static external field was cycled in the interval{
-
H
max
⁡
,
H
max
⁡
>
0
}, such that the value of H
max
⁡ is sufficient to impose a relatively large average alignment along e
x  (
e
y
) for H
=
H
max
⁡ (H
=
-
H
max
⁡), corresponding to η
~
1 (η
~
-
1). We calculated hysteresis loops at zero temperature for a fixed value of the ordering interaction and disorder strength, represented by positive constants J and D. Both mean field and lattice simulation approaches yielded qualitatively similar responses on cycling the forcing field. The hysteresis loop was either absent or can consist of single or double loops. Single loops were obtained for relatively low ratios D
/
J. For D
/
J
~
1, double loops could emerge. In this case no hysteresis is observed at small values of H. For large values of H two loops appear in a narrow window of H values for positive and negative values of H, symmetrically placed with respect to H
=
0. However, both approaches yield quantitatively different results suggesting that mean field results well predict only qualitative shape of hysteresis loops.It is of some interest to pose the question as to why the mean field theory is not more accurate. The mean field theory has been constructed, following the work of Shukla and Kharwanlang [17], following earlier work on time-dependent systems by Strogatz et al. [28]. The heuristic philosophy of this disorder-based mean field theory resembles that of standard mean field theories in the temperature domain (e.g., Curie-Weiss theory [30], or the well-known Maier-Saupe theory [31] in liquid crystals). There have of course been extensive studies of the circumstances in which thermal mean-field theories might be supposed to be reliable, both in terms of critical exponents and numerical predictions of phase transition temperatures. In general this is true for long-ranged potentials and sufficiently high dimensions, and presumably this is the case here also, but more detailed studies are required for these disorder-based mean field theories. At this stage all we can do is to state that the mean field theories are heuristically appealing and might thus be expected to give some insight into what might be expected.We note the relationship between our work and a recent paper involving some of the present authors [32]. This paper considered the finite temperature zero field behavior of the n
=
3, d
=
3 version of the model under discussion here, subject to a number of different cooling regimes. In this case too, hysteresis and history-dependent properties are found. There is surely a relation between field-dependent and temperature-dependent hysteresis, but further studies are required to elucidate it properly. In [32], the zero-field properties are intimately connected to the nature of long distance correlations. But here, at zero temperature in a field, we are only concerned with the properties of the total system. It seems likely that having a strong-field history (essentially a field-cooled system), that long-range order is maintained until the system is reoriented. However, to be sure, further size-dependent simulations are necessary, and the finite temperature hysteresis properties could in principle be different from the zero-temperature properties even at infinitesimally small temperatures. A detailed comparison with [32] is not currently possible, as the dependence of the various properties on d and n is not known at this stage.The problem posed in this paper is a pilot study of a three-dimensional liquid crystal in a random field. The addition of a third spin and a third spatial dimension presents some computational challenges. It may change the qualitative properties of the hysteresis cycle, in the same way that it does change the qualitative properties of the low temperature phase. The physical problem is that of a liquid crystal in a porous system, in which the volume of the enclosing porous material is a small proportion of the total volume but nevertheless is sufficiently strong to localize the liquid crystal. The hypothesis is that such systems are well described by random field systems of some sort. However, much detailed study is still required to validate this hypothesis, both in terms of comparison between models and comparison with experiment. The present paper merely indicates that some features of the hysteresis behaviour of systems of liquid crystals in pores seem to occur in even quite simple random field models.We note from the points above that the primary motivation of this paper is theoretical and concerns the properties of random field models. This general problem has attracted a considerable degree of interest. The rather peculiar trace-free maximally biaxial imposed field has been chosen so that one can cycle reversibly from a positive to a negative field and back, in such a way that But one might also ask the question as to whether the question as posed finds some experimental application. In fact random anisotropy models have been applied, particularly by Terentjev and coworkers [33–37], and more recently by Lopatina and Selinger [38], to nematic elastomers and related systems. For these cases, a biaxial stress field with the symmetry discussed in this paper is a natural external field, while “shape memory” is a natural concomitant of the order parameter memory discussed here.



---

*Source: 834867-2015-02-23.xml*