lilllllljllllllllflllllfljllfll j

93014302 6

 

LIBRARY
Michigan State

 

 

 

This is to certify that the
dissertation entitled

PHASE TRANSITIONS IN TWO DIMENSIONAL
DJATOMIC MOLECULAR SYSTEMS

presented by
SAN-Y | TANG

has been accepted towards fulfillment
of the requirements for

Mdegreein I HYSICS

Mb. MMa/dk/

Mafor professor

Date /0.24. 1985

 

 

 

MSU

LIBRARIES
.—_.

hut

 

 

RETURNING MATERIALS:
Place in book drop to
remove this checkout from
your record. FINES will
be charged if book is
returned after the date
stamped below.

 

 

 

 

 

PHASE TRANSITIONS IN TWO DIMENSIONAL

DIATOMIC MOLECULAR SYSTEMS

by

San-yi Tang

A DISSERTATION

Submitted to
Michigan State University
in partial fulfilment of the requirements

for the degree of

DOCTOR OF PHILOSOPHY

Department of Physics and Astronomy

1985

ABSTRACT
PHASE TRANSITIONS IN TWO DIMENSIONAL

DIATOMIC MOLECULAR SYSTEMS

by

San-yi Tang

Phase transitions in two dimensional diatomic molecular systems are
studied in this thesis by using various theoretical methods and computer
simulations.

It is shown that realistic molecular interaction can be represented
reasonably well by a spin Hamiltonian for a rigid lattice which
corresponds to fix the mass centers or the molecules. A practical
method to derive such a spin Hamiltonian is given. Under very general
conditions, this mapping leads to an aniostropic XY model.

The effect of this anisotropy on the isotropic XY model is
investigated. It the anisotropy is present, no matter how small it is,
the long-range order can exist, and as temperature rises, the system
undergoes an Ising-like order-disorder transition instead of the
Kosterlitz-Thouless transition seen in the isotropic XY model. The
physical reason is that the dominant term in the excitation energy of
vortex-antivortex pair is now proportional linearly to the separation
between the opposite vortices instead of logarithmically in the

isotropic phase.

Naturally, the next step is to release the mass centers of the
molecules.

After solving the problem of incorparating the orientational
contribution to the deformation of lattice within a constant-pressure
molecular dynamics scheme, the phase transitions in a two-dimensional
diatomic Lennard-Jones molecular system consisting of 400 molecules are
studied. The computer simulations show that the coupling between the
orientational and translational degrees of freedom drives the
orientational order-disoder transition first order. Topological defects
play an important role in the driving mechanism for this transition.
Through this transition the system transtorms trom a terroelastic phase
to a paraelastic phase. Simulations also suggest that the melting
transition of this system is first-order as has been seen in 2D

monoatomic systems .

ii

DEDICATION

To my contemporaries in China:

WHO suffer,

Who struggle,

Who conquer.

iv

ACKNOWLEDGEMENTS

I would like to thank professor S.D. Mahanti for his invaluable help
and guidance as my thesis advisor. His encouragements made this work
possible and his theoretical insight made it a pleasant expedition. The
enthusiasm he devotes to the research will be a life-long inspiration
for me. Despite the sincerity of all these gratitudes, yet only the
time coming in future could tell the true extent of his influence.

I would also like to thank Professor G. Kemeny and Professor
T.A. Kaplan for many instructive discussions and their constant interest
in my work.

I am grateful to Professor D. Bruce, Professor J. Kovacs and
Professor D. Stump for their service on my guidance committee and the
very helpful comments on this thesis.

I wish to thank Dr. R.K. Kalia, Argonne National Laboratory, for
stimulating discussions when I was visiting Argonne.

The faculty of the Department of Applied Physics, the Shanghai
University of Technology and Science, are greatly appreciated for their
hospitality when I was in China in 1984.

Finally, I would like to emphasise that without the commitments and
supports of the all members of my family, this work would never be
produced. Especially, I feel very indebted to my wife, Renqiu, and my
daughter, Yan, for those years I have had to leave them in order to

complete my education.

TABLE OF CONTENTS

LIST OF TABLES

LIST OF FI

CHAPTER 1

CHAPTER 2-

II.

III.

IV.

CHAPTER 3

II.

III.

IV.

VI.

VII.

GURES

Two Dimensional Molecular Systems:

a General Introduction

Lennard-Jones Molecules on Two Dimensional
Lattice: a Model Anisotropic XY System
Introduction

Effective Spin Hamiltonian

(Expansion Method)

Effective Parameter Procedure (EPP)
Results and Discussion

Phase Transition in

Anisotropic Planar Rotor Systems
Introduction

Spin-wave Approximation
Migdal-Kadanoff Real-space
Renormalization Group Procedure

Monte Carlo Renormalization Group
Simulations

Vortices and Strings

High Temperature Series Expansion

Summary

vi

viii

ix

10

IO

14

23

25

30

30

36

39

44

50

58

60

CHAPTER 4 Phase Transitions in

II.

III.

IV.

V.

Diatomic Molecular Monolayer

Introduction

The Ground State

The Constant-Pressure Molecular Dynamics

A) Constant-pressure Ensembles and Stress Tensor
for Multiatomic System

B) Numerical Solution of the Equations of Motion

C) Aging and Generation of
the Equilibrium Ensembles

D) Simulations of the Lennard-Jones System

The Perro-Paraelastic Phase Transition

A) Lattice Anisotropy Parameter and
Orientational Order Parameter

3) Radial Distribution Function

C) Energy and Density

D) Orientational Diffusion Coeffiecient

E) Strain Fluctuation and Elastic Constant

P) Discussions

Summary and Comments

LIST OF REFERENCES

vii

62

62

69

72

72

78

82

84

87

87

89

90

90

91

93

108

114

LIST OF TABLES

Table 2.1 Parameters of effective spin Hamiltonian

Table 2.2 Low-temperature order parameter

and transition temperature

Table 3.1 Transition temperature obtained by MKRG and MCRG

viii

27

28

49

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

2.1

2.2

2.3

2.4

3.1

3.2

3.3

3.4

3.5

3.6

3.7

4.1

4.2

LIST OF FIGURES

The relative positions of two molecules
Possible ground state configurations

Ground state configurations in parameter space
Monte Carlo results

Approximate Ising-like interaction

RG flow

Specific heat and susceptibility

MCRG results

Results of quenching studies

Possible configurations of VA pair

Density of VA pair vs. reciprocal of temperature
Ground state of 2D diatomic molecular system

Vortex on triangular lattice

ix

13

21

22

29

42

43

47

48

55

56

57

67

68

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

4.3

4.4

4.4

4.5

4.6

4.7

4.8

4.9

4.11

4.12

Lattice anisotropic parameter and

Orientational order parameter vs. temperature

(a) Radial distribution function at T=0.36

(b) Radial distribution function at T=0.38

Energy vs. temperature

Density vs.

temperature

Orientational autocorrelation function

Orientational diffusion coefficient

vs. reciprocal of temperature

(a) Configuration at T=0.36

(b) Configuration at T=0.38

Orientational potential well

Configuration of quench study

Radial distribution function at T=0.7O

Configuration at T=O.70

97

98

99

100

101

102

103

104

105

106

107

112

113

CHAPTER 1

Two-dimensional Molecular System:

a General Introduction

In condensed matter physics, the rapid growth of research on
two~dimensional systems is a notable phenomenon in the past ten
years.(1'1) Indulging in the fascinating new dimensions opened by modern
technology and science, we are reluctant to recall the historical fact
that for millions years the human being had been a poor creation
confined to a two-dimensional manifold-the surface of the earth. So it
is no wonder that an essentially two-dimensional theoretical
system--celestial mechanics, reached its maturity first in history among
all branches of the physics, and 2000 years ago the ancient Chinese
astronomers already thought that the planets were traveling in a 2D

circular orbit.(l'2)

However, at the microscopic level, the realization
of true 2D systems has been achieved only in the recent decade, though
physicists long before have known that some kinds of layered compounds
can be described quite well by 2D theories. A prominent member of the
class of 2D systems is the molecular monolayers physisorbed on graphite
substrate. The commercially available graphite materials with
high-quality surfaces have helped the advance of experiments in the
70's; these experiments have found close correspondence between the
measurements and expected 2D behavior in many cases. Extensive research
has been carried out on relatively simple systems like noble gas

(1.3)

molecules Ar, Kr, Xe on graphite, and N molecular monolayer which

2

is commensurate with the graphite substrate.(l°4) The mutual stimulation
between experiment and theory with computer simulation as a powerful
weapon, naturaly extends the frontier further forward to more

whose monolayers are incommensurate
(1.5) (1.6) on

complicated molecules such as 02

with the graphite substrate. Methane, ethane, and ethylene
graphite substrate are of great current interest. For these systems,
richer phase diagrams have been unveiled.

In this thesis, I have used a realistic interaction
model-—Lennard-Jones potential to investigate the properties of a
diatomic molecular system, which resembles closely the 5-phase of the
oxygen monolayer on graphite below 40K. In this phase, the oxygen
molecules lie flat on the substrate forming an almost centered
rectangular structure which is incommensurate with the underlying carbon
triangular lattice.

In the middle seventies, McTague and Nielson(l’7)

conducted neutron
scattering experiments to explore the phase diagram. of an oxygen
monolayer on graphite in the temperature region 4.2K<T<50K, and coverage
region 0.7<p<2.6. (The coverage 1 corresponds to a V3x¢3 monolayer
which is the structure of Kr molecular monolayer on graphite. Although
we know the oxygen monolayer never shows this structure, people still
use this unit for reasons of convenience or tradition.) At high
coverage, p>l.6, they observed two phases: above ll.9K they saw a single
diffraction peak at reciprocal lattice position Q=2.204A-l, below 11.9K
this peak split into a doublet at 2.17 and 2.3011’1. A peak at

l

1.14A- which was ascribed to magnetic ordering, was also observed.

They labelled the low-T phase as a, thought it had the same positional

and magnetic structure as the densest packed planes of bulk a-O The

2.
other phase, labelled B, was deduced to have the same triangular
structure as the basal planes of bulk 3-02. 'These two phases are

relabelled by Stephens et al.(1°9)

as e and 5, who argued that there was
no clear connection between the phases of bulk and adsorbed oxygen. The
transition between these two phases has an interesting and yet
unfinished story by itself, my interest in this thesis, however, is on
the low-coverage 6-phase. When p<l.6, McTague and Nielson found a
single peak at 2.15A—l, they thought that the 5-phase had a triangular
structure. However, it turned out that a peak at 1.57;;-1 was missing,
and was discovered‘in the x-ray experiment performed by Stephens et

1 was in fact a doublet

al..(1'9) They also found that the peak at 2.155'
with two closely separated, resolution-limited peaks. It was inferred
that the 6-phase had a centered parallelogram structure with side
lengths a=3.l7-3.25A, b=7.97-8.85A. The axes of 02 molecules are
parallel to the substrate plane. Independently, by using a pattern
search program to calculate the minimum of potential energy, Etters et
al.(l'1o) predicted the existence of this phase with a=3.32A and
b=8.34A, in very good agreement with the experiment. Later the
low-energy electron diffraction (LEED) experiment carried out by Toney

et a1.(1'll)

confirmed this discovery. Now there is a true 2D diatomic
molecular system whose phase diagram is very attractive for its
potential richness.

The two-dimensional solids are rather different objects than
ordinary three-dimensional crystals and the simple idea of long-range

density-density order is not satisfactory to characterize them. In the

thermodynamic limit, the Debye-Waller factor describing the thermal
vibrations diverges in 2D at any finite temperature. Kosterlitz and

Thouless(l'12)

proposed a different definition of long-range order for
2D solids depending on the overall properties of the system rather than
the pair correlation function. They called it ”topological long-range
order". This new criterion for solids is based on the supposition that
the 2D system has enough short-range order so that the local crystal
structure can be defined, then if there are no free dislocations
present, the system is rigid, while if the free dislocations are
present, an arbitrary small shear stress will cause the dislocations to
move and the system will have a response characteristic of a viscous
liquid. Thus the presence or absence of free dislocations determines
which phase the system has. The great advantage of this theory is that
the topological defects are much easier to deal with in 2D than in 3D

because now they are point defects. Nelson and Halperin(l°13)

developed
this idea, by defining an order parameter for bond orientation. They
argued, for spherical molecules on a smooth substrate, the melting from
solid to liquid might take place in two steps; there was an intermediate
liquid-crystal (hexatic) phase with exponential decay of translational
order but algebraic decay of six-fold orientational order. Similiar

conclusion has been reached by Young independently.(l'1‘)

This
proposition invokes a long controversial problem about the nature of 2D
melting. Contrary to the experimetal evidence and the results of

computer simulations(1'15)

for first order melting, this theory predicts
two continuous transitions for 2D melting. Later, Ostlund and

Halperin(1'16) applied this dislocation-mediated 2o melting theory to

anisotropic layers possibly consisting of rodlike molecules, and
concluded that the melting transition might also go through two stages,
but the properties of the intermediate liquid-crystal phase was more
complicated because now the system has two types of dislocations.
However, there are attempts to unify these two sides. Saito proposed
that the order of the 2D melting transition might depend on the core
energy needed to produce a topological defect pair, for large
core-energy the transition is continuous, and first-order melting is

(1°17) For 2D diatomic molecular

found in small core-energy case.
systems, the presence of the additional orientational degrees of freedom
which is absent for spherical molecules, can introduce other types of
possible new phases. Although the true translational long-range order
is absent, the long-range orientational order can exist due to the
anisotropic nature of the intermolecular interaction. Thus, between the
solid and isotropic liquid phases, a plastic phase with exponential
decay of orientational order but algebraic decay of translational order
can emerge. The properties of the 2D plastic phase and the
ferro-paraelastic phase transition are almost untouched in the
literature. In addition, a liquid-crystal phase is also possible for 2D
diatomic molecular systems under certain conditions.

To describe such a molecular system, we ought to know the
intermolecular interaction potential. The most direct way to determine
it is to solve the electronic SchrBdinger equations for all relevant
values of the molecular coordinates. However, to attain the accuracy

required for meaningful thermodynamic calculations it is necessary to

take into account the electron correlation effects, a task which has not

been performed by completely ab initio quantum-mechanical methods. The
cohesive force which binds nonpolar molecular crystals and fluids are
largely intermolecular electron correlation effects,(1°18) and are not
included, for example, in a Hartree-Fock calculation. Various attempts
have been made to construct an intermolecular interaction from additive
atom-atom potentials which are assumed to be transferrable from one
molecule to another with some success.(l'19) As mentioned above, by
using an additive atom-atom potential consisting of a Lennard-Jones part

and a quadrupole-quadrupole part, Etters et al.(l’10)

successfully
predicted the 5-phase even when the earlier neutron scattering
experiment seemed to contradict their results. Because oxygen has a
relatively small quadrupole moment and the magnetic coupling constant is
also small compared to the intermolecular dispersion

potential,(1‘20)

(though it is expected to play an important role at the
e-S phase transition in the high-coverage region), we did not include
them in our calculations. The absence of a commensurate phase shows
that the interaction between the adsorbed 02 molecules and the graphite
substrate is weak, we believe that its main role is to make the molecule
lie flat on the substrate in the 5-phase. The substrate corrugation
potential can be integrated over the surface resulting in a constant
term which can be added to the internal energy of the system. The
molecules are treated as rigid rotors. It is a good approximation to
regard the lattice vibrations as consisting of translations and
rotations of those rotors, since the intramolecular vibrations are

generally not excited at low temperatures considered here. The system

studied in this thesis consists of diatomic molecules interacting via an

intermolecular potential which is the sum of the atom-atom Lennard-Jones
potentials between the constituent atoms of different molecules.

A basic concept of the modern theory of ’critical phenomena is
”universality”; it claims that the nature of the phase transition in a
particular system is determined by the symmetry of the Hamiltonian
describing this system, the number of the components of the order
parameter and the dimensionality of the system. According to this
hypothesis, different model systems, even real physical systems, are
subject to the same ”universal class" if they lead to the same
Hamiltonian in the Ginzburg-Landau-Wilson (GLW) form on the same lattice
structure. To study the GLW Hamiltonian, many sophisticated methods of
statistical mechanics are at our disposal, including the renormalization
group method, the most fruitful procedure in the area of critical
phenomena up till now. For 2D molecular systems, it is desirable to
develop a method to derive an effective Hamiltonian which is simple
enough to be tractable but without losing the essential features of the
system. In Chapter 2, assuming that the centers of mass of the
molecules are confined to a fixed lattice, we show that under very
general conditions, the orientational motion of a 2D diatomic molecular
system can be described by an effective spin Hamiltonian in the
anisotropic XY form. The parameters of this spin Hamiltonian can be
found by the effective parameter procedure (EPP) proposed in this
chapter. We performed Monte Carlo simulations using parameters obtained
by EPP, the results are encouraging when compared with that obtained
from molecular dynamics simulation using the exact Hamiltonian. The

advantage of this effective parameter procedure over the conventional

small parameter expansion method is discussed in this chapter.

In Chapter 3, the phase diagram and the nature of the phase
transition of the effective spin Hamiltonian are studied by using
various methods, including Migdal-Kadanoff real space renormalization
group (MKRG) and Monte Carlo renormalization group (MCRG). The results
show that in this anisotropic XI model, long-range order is present at
sufficiently low temperatures and the ferro-paramagnetic phase
transition is Ising-like. Under MKRG procedure, the system iterates to
an Ising-like system, and MCRG gives the correlation length critical
exponent a value close to that of 2D Ising system. The special feature
of this model is that it supports vortex excitations while possessing
the simplest domain-wall pattern. A quench study gives the first
simulation picture which clearly shows the interaction between these two
types of excitations. We also give an energy-entropy argument to
explain why the Kosterlitz-Thouless transition observed in isotropic XY
model gives way to an Ising-like transition when the anisotropy is
present in our model.

The coupling between the different degrees of freedom can
drastically change the nature of phase tansitions. In Chapter 4, we
release the mass centers of the molecules. By minimizing the potential
energy and searching for the zero point of the internal stress tensor,
we find that the ground state of this system is a centered rectangular
lattice, and its spacings are very close to those of the 6-phase of
oxygen monolayer on graphite if we use the parameters appropriate for
oxygen molecules. We conduct constant-pressure (force/length) molecular

dynamics (MD) simulations to explore the phase diagram along the zero

pressure line. A ferro-paraelastic phase transition is discovered at
20.6K, the long-range orientational order is lost, the centered
rectangular lattice at low temperatures becomes triangular. The
discontinuities found in the orientational order parameter, the lattice
anisotropy order parameter, the energy, the density and the
orientational diffusion coefficient indicate that the transition is
first order. The physical mechanism of this first order transition is
discussed in light of the role played by topological defects. At 38.0K,
this plastic phase with triangular structure melts to a liquid. Again,
the discontinuities of the energy and the density suggest that the
transition is first order. In addition, the MD method used in our

simulations is also documented in this chapter.

10
CHAPTER 2
LENNARD-JONES MOLECULES ON A THO-DIMENSIONAL LATTICE:

A MODEL ANISOTROPIC X! SYSTEM

1. INTRODUCTION

A clear understanding of the observed structural phase transitions
in 2D molecular ”solids” requires adequate knowledge of the anisotropic
intermolecular potential. The commonly used intermolecular potential
consists of two parts, one (ves) associated with the interaction between
electrostatic multipole moments of the individual molecules and the
other (vr+vdis) is the sum of short-range steric repulsion and
polarization interactions. The latter two are usually obtained from an
atom-atom potential model (standard Lennard-Jones potential between
atoms of one molecule and those of the other). I

In contrast to vr+vdis’ the contribution to ves from the interaction
between i-th and j-th molecules can be expressed directly in terms of

. A A A .
spherical harmonics Ylm(ni) and Yl,m,(n ), fii and nj being the unit

1
vectors representing the orientations of the two molecular axes. The
maximum value of l appearing in the expansion of ves depends upon the

electric multipole moments of the molecular species under study. For

11

example, the largest contribution to the intermolecular interaction for
N2 molecules comes from the latter's electric quadrupole moment which
corresponds to l=2. The repulsion and dispersion contributions to the
intermolecular potential can in principle be expanded(2'l) as an
infinite series in powers of the parameter dad/R, where 2d is the
interatomic separation of a given molecule and R is the intermolecular
separation (see Fig. 2.1). Only when a<<l, the first few terms in the
spherical harmonic expansion (l<2) make the dominant contribution and
the resulting Hamiltonian has a simple form. But for realistic systems

(2'2) d=0.15. For these

such as 02 and N2 on a graphite _substrate
systems one has to keep many terms in the spherical harmonic expansion
and the resulting Hamiltonian takes a complicated form.

To study the adequacy of the above-mentioned expansion procedure, we
start from a 2D molecular solid for which the intermolecular potential
is of Lennard-Jones form (the corresponding Hamiltonian is denoted as
H

). construct an effective spin Hamiltonian HS and compare the

23
exact ex
exact thermodynamic properties of HS

P

p and Hex’ Properties of Hex are

available from recent molecular dynamdcs (MD) simulation

(2.3)

studies. Such a comparison will tell us about the correctness of

the procedure for obtaining HS from Hex' Our ultimate aim is to apply

P
a similar procedure to realistic 2D molecular systems like 02 and N2 on
graphite which show orientational phase transitions involving both

orientational and translational degrees of freedom.(2’2)

Although the
discussion of Sec. II and III of this chapter are quite general, for
numerical studies, we have chosen a simple system where the centers of

mass of the molecules are constrained to be on a 2D square lattice and

12

the molecular orientations are confined to the plane of the lattice.
This choice was dictated by the availability of MD simulation studies in
this system.(2'3)

The outline of this chapter is as follows. In Sec. II we carry out
a straightforward expansion of the atom-atom potential in powers of a
and analyze the thermodynamic properties of the resulting Hsp' We then

discuss the adequacy of Hs to represent Hex for a values of practical

P
interest. In Sec. III, we introduce an alternate procedure which we

refer as effective-parameter procedure (EPP) to obtain Hs and compare

P
its properties with those of Hex“ This procedure is physically
equivalent to keeping only the leading terms in the spherical harmonic
expansion, but calculating the coefficients up to all orders in a.
Strictly speaking, one arbitrarily drops terms in the Hamiltonian
involving higher 1 values even if their coefficients are not small.
Results of calculations for Hsp obtained by using EPP are given in
Sec. IV and compared with the exact results to see whether one is
justified in dropping higher-order spherical harmonic contributions. A

brief discussion of the applicability of this procedure to N and 0 on

2 2
graphite substrates follows.

 

 

 

 

 

Fig. 2.1 Two diatomic molecules with intermolecular
separation R and oriented along by and 3- directions.

(A.,Az) and (A3,A‘) are the atoms associated with the two
molecules.

14

II. EFFECTIVE SPIN HAMILTONIAN (EXPANSION METHGJ)

To obtain an effective Hamiltonian describing the intermolecular
interaction, starting from an atom-atom potential, a straightforward
method is to expand the latter in powers of the ratio a. The total
Hamiltonian describing the interaction between molecules can be
expressed in terms of the potential energy of interaction between i-th

and j-th molecules, Vij' i.e.,

H = Z V.. [2.1]
ex 6,1) 13 -
where
3 an an an an
vij vA.‘A3 + A|"A4 +VA2’A3 +vAz-A‘,’ [2.2]
and the atom-atom potentials ”Inf-A; are given (for 1 =1 and m =3) by
| .
ea _ 12 _ 6
vA.-A3 - 4e[(a/r13) (a/rl3) ]. [2.3]

In Eq. [2.3], rl3=|(Ri+dfii)-(Rj+dfij)l, is the distance between atom
A1 of the i-th molecule and the atom A3 of the j-th molecule
(see Fig. 2.1) ; R1. and R3. are the positions of the respective molecular
centers of mass. Expanding vij in powers of a and keeping terms up to

0(a‘), we find that the pair interaction energy vlj is given (in unit of

46) by

Vij=v ij +i§Lt h‘[skx cos(2wi j)+ sky

2 - 2
35.5013. .[(skx sky)cos(4wij) +25

+Jij(sixij+slysjy)

Kij[(sifi -siy sjy)cos(4¢iJ )+ (sixsjy+s iy ij)sin(4\lziJ

sin(2wij)]

kxskysin(4wij)]

15

[2.4]

)1

where wij is the angle that the vector connecting i and j sites makes

with the x-axis (see Fig. 2.1) and

s =cos(291), s =sin(29j).

k3 RY
The parameters appearing in [2.4] are given by

= (a/R)12(4+2aaa2+losata‘>-<a/R)6(4+72e2+854e‘),

hk= <a/R)12(laaa2+107sza‘)-(o/R)6(48a2+9soa‘).

= [1008(a/R)12 “120(O/R)6]a4r

Jij= [3528(a/R)12 -288(a/R)6]a‘,

Kija [6048(a/R)12 -7zo(a/R)Gja‘.

In the above, R=Rij and dad/Rij. If we consider

nearest-neighbor interaction,(2'4)

[2.5]

[2.6a]

[2.6b]

[2.6c]

[2.6d]

[2.6a]

only the

take the x-direction along a lattice

bond and sum up all the pair potentials, the terms which are linear in

cos(291), sin(291) drop out due to the symmetry of the square lattice,

and the anisotropic part of the total Hamiltonian,

Hspsaznwii Vii)

16

reduces to a simpler form given by

= z - 2
H r 41x51.x s. ) + ; J(s. >

S . +5 . S .
SP 11’ «,3; 1x 38 1y W

+ Z K(s ) [2.7]

a,” ixij-siysjy
where <ij> indicates all the nearest-neighbor pairs and D, J, K are
those given by Eqs. (6c), (6d) and (6e) with Rij equal to the
nearest-neighbor distance.' Hsp has the form of an X! model which has
both single-site and interaction anisotropies.

Before discussing the finite-T properties of Hs and comparing them

P
with results of MD simulation, we would like to discuss briefly the

nature of the ground state of Hsp for arbitrary values of anisotropy and
interaction parameters D, J and K. Using the Luttinger-Tisza

(2.5)

method, we find the ground state of Hsp for different parameter

values. The ground state is always a commensurate structure, and the

Luttinger-Tisza method works.(2'6)

We also find that, depending on the
parameter values, there are six possible ground-state configurations:
two antiferromagnetic (A171, A172), two ferromagnetic (F1, F2), and two
herringbone type (H31, 882). These are shown in Fig. 2.2. In Fig. 2.3
we give the range of parameter values for which one of the above
structures has the lowest energy.

To test the adequacy of the expansion up to 0(a4) in obtaining the
parameters of H5 , we have calculated the orientational order-disorder

P

transition temperature TC of H5 using mean-field approximation and

P
compared it with the MD results for Hex' The parameters used in the

later calculations are:

17

_ -16 °
e -l3l.3x10 ergs, 0 =3.708 A,

n =0.6438 i, and R =4.2488 5.

Using these values for the parameters in Eqs. (6c), (6d) and (6e), we
find D, J and x (expressed in unit of io'leergs) equal to 39.78, 155.38,
and 238.70, respectively.

The ground-state configuration for the above values of D, J and K is
of the type AF2 which is in agreement with the low-temperature results
of the MD simulation. For T>0, we assume that the molecules fluctuate
about the orientations they have in the ground state, i.e., we assume

that thermal averages of ij and.sj are given by

Y
1}! +157 ’7

<ij> 2 (-1) [2.8a]

and
<s. > = 0 2.8b
3y [ 1
where n is the order parameter and l. , ljy are integers labelling the

3x
lattice sites. Approximating Hsp by its mean-field value HMP' where

a z_ 2, _ lix+lf¥
HMF E [4D(six 51y) +(J+K)sixn2 ( 1) ] [2.9]

and defining six=coso and siy=sin9 (note 9 is not the orientation of

molecule from now on), we find that

T¥(e*H”’5tx)
'Tr( c"""")

 

£21,005 9 “PPM 0005(99H4NJ+K)91005 9] 49
f." ”P [-4A Gauze) +43(:+D)ncose] d9

 

[2.10]

Assuming the transition to be continuous, we expand the right-hand side

18

in powers of n and calculate the mean-field transition temperature

(2.7)

TgP in the usual way. The transition temperature is given by the

transcendental equation

 

T? g 2(J'+K) 212(flc)

[2.11]
*5 Ia‘fic)

where

12n(6) = I?’cosznOexp[-4D6cos(20)]d9 [2.12]

Numerical calculation gives the reduced transition temperature
t§F=kBT§F/e=5.3. This is about a factor of 1.6 smaller than t¥P=tE§=8.4

which was obtained from MD calculations for the full Hamiltonian

H .(2.3)

ex Since the mean-field approximation neglects effects of

fluctuations, one expects the exact transition temperature to be smaller
than tfiFL Therefore, the exact transition temperature of HSP is likely
to be considerably smaller than the exact transition temperature of Hex“
This suggests that there is a serious deficiency in representing Hex by
Hsp'

To make further checks on the validity of the above expansion

procedure used to obtain Hs we have investigated the low-temperature

P

behavior of a system described by Hsp using spin-wave approximation. At

sufficiently low temperatures, we expect the deviation from the
p by
expanding 91 about its ground-state value. Using 01:9}, in one

ground-state orientation to be small and therefore simplify Hs

sublattice, and 9j=u+03 in the other sublattice, where a; and a; are

small, at low temperatures we have

19

2 I
H = 2 4D 1- 29¢ 2 - 2 1- of-o. 2
sp i [ ( l) / ] 66°{J[ ( l J) / 1
+xti-(a;+93)2/21} [2.131
After dropping the constant term and the prime index, we find that

+
H = st const,

SP
where the spin-wave Hamiltonian Hsw is given by (for J>0, K>0, and K>J)
2
H a + - '+ - . . 014
sw (23 2x 8D)§261 (Evan .msvlaJ [2 ]

The spin-spin correlation function between one spin at the origin and
the another at'? lattice site is given by

9(r) = <cos(60-9f)> = Re<exp[i(90-0?)]> [2.15]
We can calculate g(r) following the procedure used by Wegner(2'8) for

the isotropic limit (D=K=0) and obtain

 

.3 . 2' '4 2
9(r) = exPI- 5132‘: 5m 2: Y, ) 1 [2.16]

where, cf, the energy of spin-wave with wave-vector E, is given by

at = 2(K+J-4D) [l+7(coskx+cosky)] [2 . 17]
and

7 = (K-J)/2(K+J-4D) [2.18]
In Eqs. [2.16] and [2.17] the distances are measured in units of the
nearest-neighbor distance. In the thermodynamic limit, we can replace
the summation over ‘i by an integral over the square Brillouin zone.

Following the notion of Wegner, we write

#31,
2(K".T"4D)

 

9(3) = epr- f2<f) 1 [2.19]

20

1 1 v dk sinz({-?/2)

2111 3"! dkx -Il Y |+ T(c‘5fix+c°5‘1)

 

 

f2(r) = [2.20]

If |7|<l/2, the integrand has no singularity in the first Brillouin
zone. This condition is usually satisfied in the case of diatomic
molecules physisorbed on a substrate, and is true for the parameters in
this chapter. Thus, the anisotropy of the Hamiltonian eliminates the
”infrared catastrophe” which is present in the isotropic X! model. In

the limit rec, g(?) is finite. For T+0 the order parameter is given
(209)

 

 

by
n = lim[g(?)]1/2= 1 - aT, [2.21]
7-60-
where
a fi'T 1 1' dk Fax 1 [2.22]

4(K+J-4o) (210‘ a x-.. Y |+r(cosk,.+cosi,)

We refer to a as the low-temperature order-parameter coefficient. Using

the parameter of HS obtained by expansion to 0(a‘), we find that

P
a=0.074. The value of a given by MD simulation of Hex is 0.040. Thus
comparing the values of transition temperatures tc and 3 obtained from

Hs and Hex' we conclude that the physical system represented by H5 is

P P

quite soft.

As a plausible cause of this large discrepancy between mean-field
and spin-wave results for Hsp and MD results for Hex' we have explored
the effect of further neighbors which were included in the MD
simulations. For the parameters of interest this can be ruled out. For
example, when we include 24 neighbors (in a 5x5 square lattice), the

mean-field critical temperature changes by only about 0.4%.

 

 

AJ=1

 

 

 

 

F1 F2

 

 

 

 

 

 

 

f

H81 H82

 

 

 

 

 

 

 

Fig. 2.2 Various ground-state configurations of the
molecules whose centers of mass are confined to a square

lattice. AF denotes antiferromagnetic; F denotes
ferromagnetic; HB denotes herringbone. .

22

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

PHASE DIAGRAM FOR J> O
3 r .
2 - AF1 F 1
1 1..
1 1 1 1 1 1 41
DH
.1 _
.—
-2 - F’2
F
.3 a
241 1 1 1 1 1
-5 -4 -3 -2 -1 2 3 4
Fig. 2.3 Ground-state configurations for different values

of D/J and K/J for J>0.

For J<0, turn the figure upside

down but keep K/J and D/J axes in coventional directions
positive to the right and up,respectively). HB
structure can be ground state only when J=K=0, so it does

(i.e.,

not appear in this figure.

 

23

III. EFFECTIVE-PARAMETER PROCEDURE (BPP)

Because of the above-mentioned shortcomings associated with

Hsp obtained from Hex by expanding up to 0(a4), we have to investigate
the convergence of the expansion procedure. If we expand the pair
interaction up to 0(a6), we find that for the square lattice and AFl,

(2'10) the total Hamiltonian still has the

AF2, F1, F2 configurations,
form of Eq. [2.7], but now (in units of 4e)
0 . [looa(a/R)12-lzo<a/R)6]a‘
+[72576(a/R)12-3024(a/R)6]a6 [2.23]
J = [3523<a/R>12-2aa(a/R>6]a‘
+[150528(a/R)12-4800(a/R)6]a6 [2.24]
x = [604a<a/R)12-7zo<a/R)6]a‘
+[241920(a/R)12-10080(a/R)6]a6 [2.25]
Except where a is very small (a<<l), the terms in the second brackets in
D, J and K are quite comparable to those in the first brackets. Thus, a
straightforward expansion procedure is not very encouraging. A simple
estimate shows that for a configuration where the two molecules are
oriented parallel to each other, the Taylor-series expansion of the
interaction potential energy in a is given by Ha “gene”. To obtain an
accuracy of three digits of the potential energy for a=0.152, we have to
keep terms with n=7, i.e., up to 0(al4) term. For parameters D, J and
K, the convergence is even poorer. They increase rapidly for the first
few terms and then decrease slowly. Furthermore, the Hamiltonian

becomes more complicated as it contains terms which are of higher powers

in s. ,

1x and a finite-T study of such complicated Hamiltonians is

Siy

24

practically impossible.

To avoid the difficulty associated with the expansion procedure, we
use a different method to obtain an effective spin Hamiltonian. We
assume that the interaction potential can be represented reasonably well
by the form given in Eq. [2.7]. We have therefore four quantities, v0,
D, J and K, to be determined. We determine these parameters by fitting
the four ground-state configurations AFl, AF2, F1, and F2, (two
antiferromagnetic and two ferromagnetic). We call this the
effective-parameter procedure (EPP). The values of D, J and K obtained
by this procedure are found to be considerably different from those
obtained by expanding up to either a‘ or as, and are discussed in the
next section. We should point out that the above procedure is
equivalent to writing

+ = Z 2 C Ylm(ni) Yl,m,(nj) [2.26]

Va vdis :,.¢.'lm,l'm'

and calculating the leading order coefficients C by actual

lm,l'm'
numerical integration. As we shall show in the next- section, the
coefficients obtained by the two procedures are very close (not exactly

equal) to each other.

25

IV. RESULTS AND DISCUSSION

For the values of s and a given in Sec. II, the effective-parameter
procedure gives values for D, J and K which are given in Table 2.1. The
corresponding quantities obtained by direct integration (form
Eq. [2.26]) are given inside brackets in the same table. These values
are reasonably close to each other and differ drastically from those.
obtained by the expansion procedure carried out up to 0(a4). The model
anisotropic XY system (described by Hsp) gives a reasonable value for
the mean-field transition temperature (tfif) and an excellent description
of the low-temperature behavior of the order parameter obtained from the
spin-wave approximation (see Table 2.2). It gives 0.042 for the
coefficient a (see Eq. [2.21]) and the mean-field transition temperature

t¥F is 15.3. To obtain the exact transition temperature for H5 , we use

P
a Monte Carlo (MC) procedure. Since MD calculations were performed for
a 30x30 sample, we also use a same size sample for our MC simulation.
We use periodic boundary condition and heat the system starting from the
ground-state AFZ configuration (T=0). We find that after about 1000 MC
steps per spin (MCS/S), the system reaches thermodynamic equilibrium for
almost all the temperatures, studied. We find that the transition
temperature t%s is between 8.2 to 8.7, the order parameter
(see Fig. 2.4) drops sharply from 0.5 to 0.2 in this temperature range.
In Fig. 2.4 we also give the results of the molecular dynamics
calculations of Ref. (2.3). The agreement is very good if we note that
MC calculations were performed with only 1000 MCS/S, whereas in the MD

simulations the system was allowed to run uninterruptedly for several

26

thousand time steps (10000 to 20000) near the transition region.

In summary, the anisotropic XY Hamiltonian (Eq. [2.71) with
parameters obtained either by fitting to several ordered configurations
or by numerical integration of the spherical harmonic expansion
describes the thermodynamics of the actual system interacting via
atom-atom Lennard-Jones potential reasonably well. For molecules like
N2, which possesses a relatively large electric quadrupole moment, the
quadrupole-quadrupole interaction will also contribute to the parameters

J and K of Eq. [2.7] and one should add these contributions, i.e.,

J = JLJ +JQ

1:: Km +1<Q

D = DLJ

where J KLJ and D are obtained by using the effective-parameter

LJ' LJ
procedure discussed in this chapter.

27

Table 2.1

The spin Hamiltonian parameters D, J and K of
eqn. [2.7] obtained by expansion up to 0(a4) and the
effective parameter procedure (see text). All parameters

are in units of lo-léergs.

 

 

 

Effective-parameter Expansin
parameters procedure to 0(a )
D 222 (221) 40
J 524 (483) 155

K 755 (741) 239-

 

 

28

Table 2.2
The low-temperature order parameter coefficient a obtained
by spin wave approximation and the reduced transition
temperature tc for different Hamiltonians. tC has been
obtained using mean-field approximation, Monte Carlo, and

molecular dynamics calculations.

 

 

 

Model a 1:"? t'Em
. 4 a
ExpanSlon 0(a ) 0.074 5.3 x
Effective-parameter a c
procedure 0.042 15.3 8.510.2
Exact Hamiltonian 0 .040b x 8.3:t0.lb

 

 

a: Spin wave theory.
b: Molecular Dynamics.

c: Monte Carlo.

29

 

 

 

 

 

(1J3 1 L l I 1 L l l 1
1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
t
Fig. 2.4 Temperature dependence of the order parameter n.

x denotes molecular dynamics results of H , 0 denotes
Monte Carlo results of H , (obtained by usi effective
parameter procedure agscussed in Sec. III of this
chapter); I denotes spin wave results of Hsp.

30
CHAPTER 3
PHASE TRANSITION Ill

ANISOTROPIC PLAKAR ROTOR SYSTEM

1. INTRODUCTION

In recent years there has been considerable interest(3'l) in

2-dimensional (2D) or quasi-2D physical systems consisting of molecules

adsorbed on a surface. Under certain physical conditions,(3'2) the

molecular interaction vij can be represented reasonably well in the
form(3.3)

vij = v0 -Jij cos(Oi-9j) -Kij

-Aij[cos(201+zwij)+cos(20j+2wij)] [3.1]

cos(oi+9j+zwij)

where oi is related to the orientation of the i-th molecular axis and
wij is an angle describing the orientation of the intermolecular bond
with respect to a preferred anisotropy axis. A useful method to derive
such a spin representation for a realistic interaction was discussed in
chapter 2. This mapping enables us to use various theoretical methods
developed for attacking magnetism problems to the study of phase

transition in 2D molecular systems.

The simplest Hamiltonian arising from vij is of the form

31

= - 9.- . - .+ .+ .. .
H J‘Ei)cos( l 93) Kai'ipcosw‘1 0] 4’13) [3 21

The above Hamiltonian was used to study the properties of H2 molecules

on a graphite substrate for particular values of J and K.(3'4)

However,
a detailed study of the properties of H given in [3.2] in the entire
parameter space is lacking and should be extremely interesting because
of the growing attention to the molecular overlayer systems physisorbed
on solid substrates. It is our purpose to present such a study in this
chapter. Furthermore, this Hamiltonian by itself is well worth
scrutinizing theoretically from the point of view of the role played by
topological excitations in the phase transitions of 2D systems. Great
effort was given to explore the driving mechanism of the phase
transition in the present model.

For simplicity, we choose J>0, K>0, e=K/J<1 and assume the rotors to
be on a square lattice. On a square lattice, the angle ¢ij is a

multiple of 21, so it can be dropped and H is now written as

H = -J 2 cos 9.-0. -K 2 cos 0.+0. 3.3
M) (l 3) «w <1 1) [ 1

In this form, for K>J case, by replacing oi with -91 in one sublattice,
we can interchange the roles of J and K in [3.3] and map it to J>K case.

The Hamiltonian given in [3.3] reduces to the well-known isotropic
planar rotor (classical HY model) model for H80. The 2D xx model has a
very instructive history. For 2D isotropic Heisenberg model, which
covers the 2D XY model, Mermin and Wagner(3'5) have rigorously proved
that the long-range order can not exist at any finite temperature. But

(3.6)

Stanley and Kaplan pointed out the possibility of the existence of

a low—temperature (low-T) in which, though LED is lost, the spin-spin

32

correlation decays so slowly that the susceptibility becomes infinite
(and remains so throughout this low-T phase), and gave evidence of a
finite temperature phase transition for the 2D Heisenberg model, based
on the analysis of high temperature series expansion for the
susceptibility. Stanley gave similar evidence of such a transition for

the X! model. Hosterlitz and Thouless<3°7)

clarified the nature of the
low-T phase by introducing the concept of ”topological order”. They
contended that the low-T phase is characterized by the power-law decay

of the spin-spin correlation, as proved by Wagner(3'8)

with spin-wave
approximation. However, the presence of the ”topological defects”,
i.e., tightly bound vortex-antivortex pairs, modifies the harmonic
picture. Now the local minimum, from which the small angle deviations
excited by spin-wave are measured, not only has some kind of short-range
order, but also contains vortices which interact with each other via a
2D coulomb potential logarithmic with respect to their separations. At
the transition the pairs unbind and it leads to a new phase in which the
spin-spin correlation decays exponentially. Bearing this new idea in

(3.9)

mind, Jose et al. have studied the effect of symmetry-breaking

field terms of the form hpcos(p0) on the isotropic xx model (JKKN
model). They found for p>4, the system has two phase transitions; the
one seen at lower temperature is conventional and the other at high
temperature is KT-type. Between the low-T p-state ordered phase and
high-T paramagnetic phase, there is an intermediate XY-like phase. If
p=4, the system only has one phase transition but with non-universal

critical exponents. For p=2,3, the system shows one phase transition

but its nature is unclear when hp is small. Similar picture is found in

33

the isotropic p-state clock model (zp model). For 2p model, Einhorn et
al.(3'10) have argued that in the low temperature ordered phase the
vortex-antivortex (VA) pairs are bound to each other with strings which
are domain boundaries constructed by the links in the dual lattice. In
the original lattice, the orientation of spin changes by almost the same
amount when one crosses any link of a particular string. The strings

become floppy at certain temperature T this results in the loss of LRD

1.
but the VA pairs are still bound to each other through a logarithmic
potential. It leads to a KT-like phase which is then destroyed by the
unbinding of VA pairs at a higher temperature T2.

The theoretical advance in this area has been helped by the progress
in computer simulations. Due to the lack of direct experimental
verification, the computer simulation is an indispensable source of
information. For example, the beautiful Monte Carlo results of

Tobochnik and Chester(3'11)

give a convincing proof of the HT theory.

Motivated by the above development, the immediate questions we would
ask when faced with our model described by [3.3] are a) is there any
phase transition in the system, if yes, then, how many; b) what is the
nature of the phase transition, if the system could at least have one.
Especially when the anisotropy is very small, does it have a
significient effect on the KT transitionz The quest for answers to the
above questions concerns the rest of this chapter.

A special feature of our model described by [3.3] is that it allows
for vortex excitation and yet has the simplest, namely 2, domain pattern

due to the 2-fold degeneracy of the ground-state. In particular, we

investigate a) the effect of the interaction anisotropy (K¢0) on the

34

Kosterlitz-Thouless (KT) transition;(3'7)

b) the physical picture
underlying the order-disorder transition for 2D anisotropic rotors; and
c) the interplay of vortex and domain-wall excitation. The study of
this interplay as a function of anisotropy is very important because of
the role the domains show in the two phase transitions seen in the
2P model for p>pc>4. However, for the 2p model, the 2-domain system
(982) does not have vortex excitations, and when vortices do appear
(p>4), the domain-wall structure is quite complicated for simulation

studies.(3'12)

Therefore the 2p model is not ideally suited for
studying the interaction between domains and vortices.

The remainder of this chapter is organized into 6 sections. In
Sec. II we give the low-temperature spin-wave results to show the
existence of the long-range order. Sec. III is devoted to the
Migdal-Kadanoff real-space renormalization group (MKRG) procedure. In
the low-temperature spin-wave approximation, our model is similar to the

(3‘9) when the

JKKN model with p=2. However, for the p=2 case,
single-site anisotropy is extremely small, i.e., the most interesting
situation, Jose et al. did not give a conclusive result regarding the
nature of the phase transition--whether it is KT type or something else,
within the framework of Migdal-Kadanoff renormalization group. By using
our real-space MKRG, we were able to show, even when H is as small as
H/J=0.0001, our model iterated to an Ising~like system so the phase
transition should be Ising-like. To substantiate our MKRG results, we
performed Monte Carlo renormalization group (MCRG) simulations and the

results are presented in Sec. IV. MCRG shows that the phase transition

is Ising-like and the critical temperatures given by MKRG are in

35

excellent agreement with that obtained in MCRG. In Sec. V, applying a
quench process, we are able to give the first computer simulation
picture which shows the interaction between the vortices and domain-wall
excitations. Unlike what one may intuitively imagine that all the
strings from one vortex terminate at an opposite vortex thus completing
a VA pair, we find that the strings from one vortex tend to bifurcate in
our model and to meet more than one vortex with opposite sign. Then we
provide an energy-entropy argument to explain this phenomenon and why,
as long as the anisotropy is present, the phase transition shifts from
RT type seen in an isotropic XY model to Ising-like in our model.
Sec.VI contains an additional check for our results by employing high
temperature series expansion approach. Finally, a short summary is

presented in Sec. VII.

36

II. LW-TBIPERATURB SPIN-WAVE APPROXIMATION

The ground state of our model is 2-fold degenerate, the spins orient
either in the positive x-direction or in the negtive x-direction. At
sufficiently low temperature, we can assume the spins only deviate a
small angle from their zero-temperature positions. We start from one of
the ground-state configurations in which the spins orient along the
x-direction, and denote the angular deviations as 0'. We rewrite [3.3]
as

a = -J z [l-(of-91)2/2] -x r [l-(9f+ei)2/2] [3.4]
«,9 l 3 (6.2) l 3

(3.8)

Following Wegner, after a Fourier transform, we find the spin-wave

excitation spectum is given by

. = 2(J+K) [ 1- 1““

Gk m (COSkx‘PCOSkY) 1 [3 . 5]

If Kao, then (J-K)/2(J+K) will be always less than half, thus the
excitation energy will never be zero if the anisotropy is present. It
means that the ground-state structure where we start from is stable uder
spin-wave purturbation. The system can possess LRO at sufficiently low
temperatures. In contrast, there is no energy gap between the
ground-state and excited states for the isotropic XY model and long
wave-length spin-wave excitations destroy LRO at any finite temperature.

The order parameter n=<cosoi> can be calculated from the correlation

function

37

RQT Z I ‘C05(;‘?)
0’ i (5‘

 

9(?)=exp[-

In the thermodynamic limit we can replace the sum by an integral over

the first Brillouin zone,

 

-‘ - _*al’_ .
9(r) exp[ 4”“) fit) 1 [3.6]
where
a v _ .‘
f(r)= 1 II ' “S“ r) dk dk [3.7]

(21111 -1. prawns-cost» 8 Y

and 7=(J-K)/2(J+K). If K¢0, 7 will be less than l/2, so there is no
singularity in the integrand of [3.7]; thus the anisotropy eliminates
the ”infrared catastrophe” seen in the isotropic XY model. In the
thermodynamic limit, the relation between the order parameter and the
correlation function is
2 . a
n = 1m 90') [3.8]
r906
For very large r, we can drop the fast oscillating term cos(k-r), then

by transforming [3.7] into an elliptic integral, we were able to derive

an expression for the order parameter

 

- _A_.1_;_~ (271-1121 22n
"’3“ 41 (1+I€)2nu (21m: 1“ }
, in 1 1 2 9 4 6
-exp 4]. (ME)! [1+ 4a+ 64a +O(a )1}. [3.9]

where

'38

I-IE )2
[+12 °

 

e=K/J, a=(

The convergence is satisfactory for H not too small, for example, if K
is 0.1J, the expansion parameter a2 has a value 0.0728. When e=0, i.e.

a=l, the series in [3.9] will diverge. In fact, in the small K limit,

 

the leading term in this series is lne.

at
[3.9] was used to check the Monte Carlo results at low temperature;

the agreement is good below a third of the critical temperature.

39

III. HIGDAL-RADANOFP REAL-SPACE RENORMALIZATION GROUP PROCEDURE

To get a general idea about the phase diagram for H given in [3.3]

"e have used the Migdal-Kadanoff(3'13)

real-space renormalization group
procedure. We first move the vertical bonds and decimate the horizontal
ones by integrating the partition function directly, and then move the
bonds horizontally while decimating the vertical bonds. After these two
operations, we take an average to restore the symmetry of the
renormalized Hamiltonian. The above procedure defines an iteration in
our calculation. We use a 60x60 matrix to store the interaction
V(Oi,9j), after the initial values being input, all the numerical

evolution can be done on it, the renormalized interaction can be

obtained by
' I
V (01,02) = -1n [1.d03exp[-V(al,03)-V(02,03)], , [3.10]

where the absorption of the temperature into the potential is
understood.

Even for a Hamiltonian initially without single-site field terms of
the form hpcos(p0), if Kao, such terms will be generated in MKRG
procedure. This is a reason why we can drop field term in [3.1] to get
[3.2] without sacrificing the generality of our model. The treatment of
this on-site field terms in MKRG procedure is non-trivial. We did not
move them along with the bonds in view of the fact that when the

Hamiltonian consists of on-site field terms only, it does not change

after decimation. At high temperature our Hamiltonian does iterate to a

40

form where only the on-site field terms survive. So, after the

iteration, we calculate the Fourier coefficients
1
Cm= I; doldezcos(m91)V(ol,02)

1
Sn- I! 69 d9

-1 l 25in(n91)V(91,62)

up to m,n =6, which determine the strength of the on-site field terms,
and keep them in the sites when we move the bonds which is now treated
as only associated with the interaction. Bxplicitly, now eqn. [3.10]

can be rewritten as

v'<ol.02>+ ghgiqwlmejn [3.11]
= -ln 1: d03exp[-V(61.93)-V(92.93)- r h,q(o3)]

. 2 h [ (a 1+ (9 )1.
9 4 q 1 q 2

where q(0) is of the forms cos(p9) or sin(p0). We look upon the on-site
field terms generated in the MKRG procedure as some kind of average
field produced by the decimated neighboring spins. We believe that this
treatment of the on-site field terms improves the earlier calculations
which effectively moved both the on-site field term and the interaction
term.

As a check of the reliability of the above procedure, we have
applied it to the JKKN model. When hp=0, i.e., the isotropic XY case,

there are two regions of J, in one the coupling J iterates to zero

directly, but in the other J first increases to a certain value, then

‘1

decreases very slowly to zero. It indecates that during the iterations,
certain other kinds of interactions with same global symmetry as J are
produced though they are very small. We estimate that the turnout point
is J=l.l, thus Tc=0.91, reasonable in comparison with the values
obtained by other methods. For example, Tobochnik and Chester derived

For p=6, it appears there are three regions, one in
which J goes to zero but h6 has a fixed value, in the second both J and
h6 grow to infinity, and in the third (intermediate region), J first
increases and then decreases slowly to zero. This is consistent with

the findings of Jose et al.(3°9)

discussed in Sec.I Introduction.
For our model, we find that as long as K¢0 (K/J20.0001), only after
10 to 20 iterations, the renormalization Hamiltonian H? can be

represented by a simple form (see Fig. 3.1)

H = -Jx Z cosoicosej -A Z cos(291) » [3.12]

with a small correction of the form ZApcos(poi). p>2. There exists a
temperature Tc such that for T<Tc, both J' and A' iterate to infinity
(T=0 fixed point); for T>Tc, A' iterates to a fixed value while J'
approaches zero implying that the system iterates to a non-interacting
Ising spin system. The form [3.12] is Ising-like, it suggests that the
nature of phase transition could be Ising-like. The transiton
temperatures obtained from MKRG are very close to those given by Monte
Carlo renormalization group simulations (see Table 3.1). The projection

of the RG flow on J-K plane is given in Fig. 3.2.

 

 

 

 

 

 

 

0.4

I
~
on: ,.

 

 

 

Fig. 3.1 Comparison of the approximate interaction (solid
line) given in eqn. [3.12] with the renormalized
interaction after 10 iterations. One spin is fixed at
080, the abscissa is the orientation of another spin.
Top: the initial parameters are Jsl.l, K80.0001J (iterate
to infinity). Bottom: the initial parameters are J=l.2,
K=0.0001J (iterate to zero).

43

 

 

 

 

- 4.0
“008
~26
-a4
\
\
\
\
\
\\
\ ~01
\
\
\\\
I I l [ ‘~‘~'-
0.2 0.4 0.6 0.8 10
T
Fig. 3.2 The projection of RC flow generated in

Migdal-Kadanoff real space renormalization group
procedure on the J-K parameter plane. The dashed line
indicates the boundary of two regions, in one region the
flows go to the zero temperature fixed point, in the
other they go to the infinite temperature fixed point.

44

IV. m CARLO RENOIU‘IALIZATIOR WP SIMULATIONS

To study the thermodynamic properties of H given in [3.3] and
support our MKRG results of the previous section, we have performed
Monte Carlo simulation of NxN systems with N=l6 and 32. To prepare our
samples, we start from a random configuration at the reduced temperature
T*=T/kBJ=2.0, which is supposed to be well above the transition
temperature, where kB is the Boltzmann constant. From now, we will
always talk of the reduced temperature and drop the * mark except in

Sec. VI. Using the standard Metropolis(3'14)

procedure, we let the
systems run at this temperature for a few thousand Monte Carlo steps per
spin (MCS/S); we then reduced -the temperature gradually with a
temperature step less than 0.2, until T=0.2 was reached and the order
was well developed. For the small K case, the systems were cooled to a
temperature equal to K. After the cooling process, we heated the
systems back. At every temperature stop, the number of MCS/S were
usually 5000-8000; however, near the the critical temperature,
l.lxlo‘ MCS/S were discarded and 1.23110“ MCS/S were used to compute the
thermal averages.

Calculations of the thermodynamic quantities such as average
magnetization (n), specific heat (C) and susceptibility (x) indicate
that the system shows only one phase transition for the three values of
K/J chosen in our simulations (0.01, 0.1, 1.0). We have found that C
and x increase rapidly near the critical temperature, furthermore they

tend to peak at the same temperature (see Fig. 3.3) suggesting that the

transition is Ising-like. In contrast, for K=0, it is

45

ed(3.11)

believ that C peaks at a temperature slightly higher than

Tc=TKT where x diverges.
Since it is difficult to locate the transition temperature from the
T-dependence of thermodynamic quantities in simulation studies in finite
size systems, we have used the MCRG procedure proposed by Schenker and

To]: hcik(3‘15)

to find Tc. The essential feature of this procedure is
the following: one starts from two systems of different sizes, say 1024
spins and 256 spins, the block spin is obtained by summing the spins
around a plaquette vectorially, and then normalizing the sum. The
plaquettes are selected in a way so that every spin is used and only
used once to construct the block spin. One repeats this procedure until
reached the limit imposed by the size of the system, such as 2x2. Our
experience shows that the thermodynamic average on a 2x2 block spin
lattice usually has bizarre behavior therefore we didn't use them.
Next, thermodynamic quantities are calculated and matched for two block
spin lattices of the same size but originating from different spin
systems. Denoting T1 as the temperature of the originally larger spin
system and T2 as that of the smaller one; if the thermodynamic
quantities of the two block spin lattices match, we expected that the
correlation lengthes are same in these two block spin lattices, then if
AT(=T2-T1)=0, the correlation lengthes are infinite in the original spin
systems because after different numbers of iteration they are equal in
the block spin lattices, thus in this case we have a critical point; and
AT<0 is an ordered phase, AT>0 is a disordered phase. The correlation

length critical exponent u can be derived from the formula(3‘ls)

46

v= ln2/ln(l+ é—AI

d11|T1=To) [3.13]

In Fig. 3.4 we plot the nearest-neighbor correlation functions for the
8x8 block spin lattices obtained from the two systems, the larger one
(32x32) denoted by 0 and the smaller (16x16) denoted by O. The crossing
of the two curves gives the transition temperature TC and it only shows
one phase transition. The 4x4 block spin lattices give similiar
results.

In Table 3.1, we give the value of Tc measured in units of J(l+e)
for three different 6 values. The MKRG and MCRG give almost same
critical temperatures thus these two procedures mutually corroborate
their applicability for investigating this model. All three values of u
are much close to 1, the Ising value, rather than to a, the HT value.
Therefore, the results are consistent with the conjecture that a
non-zero anisotropy (Rec) makes the transition Ising-like.

However, lacking a comparison with the results for the isotropic XY
model, we are unable to say whether this MCRG procedure is sensitive
enough to distinguish the difference between Ising-like and KT-like
transitions when the sizes of the systems allow us only to perform a few
iterations, especially if the anisotropy is very small. In present
situation, the fact that val is suggestive only when it is related to
other evidences for an Ising-like transition in the model studied here.
We will return to this concern about the finite size of the systems used
in the above MCRG procedure after the calculation of the excitation

energy for vortex in our model in Sec. v.

Fig.

47

 

 

 

 

 

LnCX) C
n
o
1 L‘
41
'3
1
[1
3.
~2
2.
(1
-1
11
1
o 0
O
1
o r . . e
no as to 4 L5 '1‘
3.3 Specific heat (denoted by O) and susceptibility
(denoted by 0) obtained in the Monte Carlo simulations
for =0.1J case. Note C and x diverge at same

temperature.

.meoum>u Aunxunv

48

0C6 Awpxw—v 50...; UC—uL6uu 00.50qu 0L0 mmovuam_ C—Qm xOO—n uxm ..0.0 OED v.0 .O.—u1\¥

 

 

 

 

to» m.n>_oce one: ac.m:_«h Co co.«ucau e no co.«m_0ctoo concu.oc-~mmcmmz v.n .u_m
L. L. L.
n... O... m.O m; 0.9 0.0 m.—. 0; m.0
e 7 e _ _ _ _ _ _
gooemmxmfooux xooemxmjonx gooemmxm._.ux
. O
O O O O l
O
o oo
0
o o o
L
O
O
C
C C
. l
O
O O O
o . . .
o o o e o o e e
o O o . . . C .
ooomm oooooommm

 

 

m6

O.FN

m...

O.N

49

Table 3.1

The critical temperatures obtained by MKRG and MCRG for
three K/J values 0.01, 0.10 and 1.00, where T*=kBT /J(l+s).

The correlation length critical exponents v calculated by MCRG

 

 

 

are also presented, these values are close to "Isinggl'o' For
.9
XY model, vxy=-.(3 )
MKRG MCRG
e T' T* v
0.01 1.02 1.04 1.00:0.04
0.10 1.18 1.20 0.9710.09

1.00 1.32 1.34 1.04:0.08

 

 

V. VORTICES AND STRINGS

If we examine the configuration generated in our MC simulations, we
can see that after the phase transition takes place, there are domains
in which the spin orientations are close to one of the two ground-state
structures, and the VA pairs tend to reside inside the walls.
(see Fig. 3.5). The vortex and antivortex are defined as following: we
calculate the sum of differences of the spin orientation (in the range
from -w to a) traveling around a plaquette counterclockwise; if the sum
is 2:, we define a vortex in the center of this plaqutte, if the sum is
-22, then we have an antivortex.

To understand clearly what types of excitations destroy the
long-range order, we have made a series of MC quench studies. At a
temperature above the critical point, the density of vortex-antivortex
pairs is considerately high. However, most of them only have a few

MCS/S ”life-time”. From a topological point of view,(3°16)

a tightly
bound VA pair is essentially different from an isolated vortex. The
former can annihilate itself by rearranging spins locally, but the
effect of the latter in the system extends to everywhere no matter how
far it is away from the center of the vortex. Therefore, to eliminate
an isolated vortex, or a well-developed defect, all the spins in the
system have to be reoriented, consequently these defects would have a
much longer MC "life-time” in a quench process. After suppressing the
short ”life-time" noise, an adequate quench picture is expected to

reveal those well-developed defects.

In Fig. 3.5(a-d), we give our quench results starting from two

51

initial temperatures T>Tc and T<Tc‘ In both cases, the system was
quenched instantly to T=0.l and then monitored up to 2000 MCS/S. For
the system quenched from T<Tc, the VA pairs annihilated each other after
50 MCS/S, and after a few hundred MCS/S, the order parameter goes back
to the value corresponding to T=0.l. In contrast, the quench from the
high T (T>Tc) phase shows drastically different behavior. In a few
MCS/S, again most of the closely spaced VA pair which are trapped inside
wide domain walls annihilated each other but after a long "time” one
finds that relatively long lived defects remain in the system. They
consist of VA pairs connected by relatively sharp domain walls (strings

on the dual lattice(3'lo)

). The total magnetization is near zero after
2000 MCS/S after quench. This picture indicates that the LRD is
destroyed by the formation of domains like what heppens in the 2D Ising

model.(3'l7)

Thus our MKRG and MCRG results that the phase transition is
Ising-like are further supported by our quench study.

It can be easily shown that the energy of a single vortex in a
system described by [3.31 is same as in the isotropic xx model, i.e.,
the anisotropic parameter K does not enter the total energy of this
configuration. But the ground state energy now is -ZN(J+K) instead of
-2NJ, where N is the number of spins in the system, so the excition
energy of a single vortex becomes

a z rJlnN +2101, [3.14]
whereas in the isotropic XY model the excitation energy only contains
the first term in [3.14]. Thus in the thermodynamic limit the dominant

term in vortex excitation energy is no longer logarithmic with respect

to the size as for the isotropic XY model but linear to it as long as

52

the anisotropy is present, xto, no matter how small it is.

For finite systems, we anticipate that in order to find Ising-like
behavior, the system size should be larger than a characteristic length
for our model when the anisotropy is small. In the length scale less
than this characteristic length, the first term in [3.14] is dominant
and the system showes an XY-like behavior; however, beyond this
characteristic length scale, the second one in [3.14] becomes leading
term and the Ising-like behavior prevails. By making these two terms
equal, we can give an estimate for this characteristic length L. If
K=0.lJ, L is about 8 lattice spacing; for the K=0.01J case, L increases
to 33 lattice spacing. Remembering the largest system we used is 32x32,
we stopped at this case K-0.01J in our MCRG simulations. Although the
system may show XY-like behavior at intermediate length scale for small
anisotropy case, the nature of the phase transition is determined by the
large-scale behavior of the system. The results of our MCRG studies are
probably correct on the ground that it seems possible to iterate the
xr-like behavior (if it is present in our systems) out and to reach the
Ising-like region in our finite systems for the three ratios of H/J we
chose. The existence of a characteristic length also demonstrates
itself in our MKRG procedure. If we start with J very close to 1.1 (the
turnpoint in J axis) and very small ratio of K/J, such as 0.0001, the RC
flow will linger in the neighborhood of this turnpoint during first
several iterations, then go to K=J line in a faster pace. If H is not
so small, the RG flow goes to K=J line in a few iterations. Note if
KaJ, the second term in [3.14] is larger than the first term for any

integer N.

53

The particular type of long-lived defect structure that we see can
be understood from a simple energy-entropy argument. The excitation
energy can be reduced to be finite through managing a VA pair separated
by finite distance 1. Assuming the difference of the spin orientation
between the two domains separated by a wall is equally divided among the
strings of this wall because this arrangement reduces the excitation
energy, we find that the energy cost per.unit length (lattice spacing)

A! for a wall of width w joining a VA pair to be roughly
_ Zfl’
AE- wJ( l-cos—ay) +2Kw [3.15]

If K80, AB decreases continuously as weo, so there is no clear domain
wall pattern in the isotropic X! model; but if Kao, the minimum of AB
shifts to a finite w; for such a wall, the energy cost is proportional
linearly to the separation between vortex and antivortex instead of
logarithmically as for the isotropic case. At finite temperatures, such
an energy cost can be balanced only by an entropy term which is also
proportional linearly to the size in order to minimize the free energy.
This entropy could not come from free vortices, because in 2D vortex is
a point defect whose entropy is proportional logarithmically to the
.size. So a mechanism different from the unbinding of VA pair ought to
be introduced. The entropy associated with the domain wall of length l

is linear in 1, therefore to minimize the free energy, it becomes
favorable to produce these domain-walls at a certain critical
temperature and their floppiness, i.e., the tension of the string goes
to zero, is the cause of the phase transition. It is also the reason

why the phase transition becomes Ising-like in the presence of the

54

anisotropy K term.

An interesting feature of Fig. 3.5 is that the strings, contrary to
what one may intuitively imagine, do not start from one vortex and then
terminate at an opposite vortex so to complete a VA pair. Rather they
like to connect as many different vortices of opposite sign as possible.
Thus the opposite vortices can align alternately along the strings to
build the walls. We show two possible domain-wall configurations
connecting a VA pair in Figs. 3.6(a) and 3.6(b). The first 3.6(a) is
discussed above, in the second 4.6(b), the walls become thinner but the
total energy cost is same as in 3.6(a). Since the total wall length is
now doubled and the second configuration has more entropy thus it is
expected to have a higher probability to be formed at sufficiently high
temperatures, but in this configuration the strings also have a higher
probability to meet more vortices than one, and it is indeed the case
seen in our simulation.

To comprehend why a KT-like phase transition is absent in our model
after the Ising-like phase transition takes place, we plot the log of
the vortex pair density vs. the reciprocal of temperature.
(see Fig. 3.7) The slope which gives a measurement of the chemical
potential needed to produce a VA pair reduces drastically at the
critical temperature. This is similar to what happens at the HT

(3'11) in Fig. 3.7 we

transition temperature in the isotropic XY model,
insert the result of the Monte Carlo simulation performed by Tobochnik
and Chester for isotropic XY model. We believe that above the
Ising-like phase transition the effective VA interaction is too weak to

allow an intermediate XY-like phase as in zp(p>4) case.

 

 

«fl

     

  

 

  
  

 

 

 

 

 

 

_(a_)

 

TSUO

 

 

 

,P'P'V'PPIYVYPI'VVP‘PVPQRVQ1'V'I m 9 ‘V‘V'I‘PV'Q'YI
7"P'73'I'V"90..""9'"...Or"??? ( ‘ 'Y'O‘If‘IVIIVO
I'73'Q99"PYYQQIPIIQ'VP'VR'V'VQ' I o.""""'""
7'77'9999'9799797'78""""!"' 0. '77'17797'79‘I
'9'97972'791P'91999970'74'IQQ.or 9. 7'199'79772819
9,7,97'790717'977'V'978899'!!!9.? 1. YQQIQVPPVVVQQV
I'QQVPQ'YPVQV'V'IVQ'V‘VVQIOPIQVP Q PVVVPPDDle'ee!
9'9'7‘3977990’PVPPOQVVQPVPPYQV'I a. '77,?! .09
'PI'I'R'P'P'"""""VYIQV‘VVY' 24" ‘O o 4441-
"IPI'Q'VVVIV'VY‘VVY'QIQYVVC.IVYV 4 or? 45444444
V'VQ'V'VVVIP'999177'999179!1"oV 4 '7 44444444
9'299297979091797797VIVVVVC1'1.-V 4 V! 4I444444
97"9'07'0'VVQVVV'VVQQVVVQIVIVQ- 4 9.. 44444404
P'V'VVPYDVYYVQ'7'9"1\V7'V'oPO." 4 7' 44644444
"'.I""'."I"I"""O'..""“' ‘ 'o “.“‘.O
"""V""9'7"“"‘0"","VI' 4 I 444444|4
""’,'V""QII'I\QP‘RQVVPVPUQY' 4 .4444444
.,"""""'.‘.“"'.\'1""""' . .“..“‘
‘...“"""‘I"'C"".""""I ‘ ...“.I.
"'.""..""‘..,""~“"""'. . .‘.a‘.“-
’1'","'"”""~N’,."""""b . I..‘.“‘.I
""""""""'""""",,"' . .Q'UO I
[
9"‘9'P'IPY‘OVVOV'P'O.IYV"'Q"IQ 4 I a . 'Veeve.
VV‘Q'VV'V'VVVVVVVVIVVIV"VVI'VQV 4 0. PO. I ‘1
7‘9997977997797V!!i!1!'7‘fiivV'VI h I or? P 1'
Q'Y'V'VVQPVQVVYP71'1977‘11‘18977 - V I! P V-
"7'?"'O."""VI""VO7'.1"..I ""' "' ‘RP'P'Y
"""""'O."'P"."“"""'II ""' "'I"""'
IV'VI'I'OVVYIVVYVVVVY'VPIV'QRI'V 8."?! '9'!!V'\VYV
uIV'VVV'!QVVVVOVDV'VVVVIVOVVIVQP I'VVV 'OIVVI'\VVC
cc'V'P'IP'V'V'OYVD'I'V'VI'VO'V'o "QV‘ V‘RQ‘VR‘V9vlev.
'VY'I'P'II.\?-'v""""7f¢"‘9'7 "Vl' '99-‘7‘\\44\\1v
4]
\Y'!o"'f§ret‘\\o."9.l"lol""fa d4.6440\'\\gfifi\"'fl'ta40“}..4UI4Q'QA
‘o‘\""’\“"\Dfi!‘uQ'ACPIA19.79 (.an'}"44s\"I.A""IIPQ"V.4.RIlVP-
..I'or9.!'\V"e\"9e\A!17F1P1PQ-‘ {FOPfiP444O‘Al'V’IJ'R-oIVQQIt'4’I-
F"!...a‘KV'PP‘JIDQQI'Y'PQsR'eVPY! A§\P"444."DQ'P\I4IO.00(‘4‘7757“!
P'P‘an!“e‘)at9p7‘fi\‘§f‘\ahoQ'Po \‘oo'!.-e'.Y«A-&"e.‘4osOQA'P‘Qavva
\l'Q’neP‘VeI'V‘PQ'Fofof‘e'I-Q'QVne \‘aa4o-eOlr'n-Oo""ow.‘lf""'¢!'!
"|‘-.o.lfio\'\11“'a‘Q'O‘et'.a‘0'. P'.c‘eb4b406n\.)'\'1."v040q774.-an.
7“:a\8\¢""Qe'a‘DI‘V'IRPQp'Ine ‘-o4.~h\u44-hao‘4\flrr.erena‘OaOe'ee
P'IPPR‘KPQIDPOQ‘R'P'PPQR‘ICAV'QP .¢.¢.'a4'.'i4440fc“§t'4."t'.4‘4500.4t
'Vfif"'n'\'P'-...c‘9P‘VPPPPcoa'a'I 44Ol‘DIP4QO-o\...4.§'\\'.qflsnu.o‘444444
VVV'QeH-OFVPP‘Qa.’.‘\!"a‘4'!A)I'V!p 411144....0'Q'Ora4.!'Qp'4pde‘54'.44un
"P‘QP'.""\1QI"A\'DAOP\I.RQ"VV 4444lidov.4(‘4!,“"n4"m4l¢e4‘4
))Q§r£9331.‘P‘K'Q!"QD'Y.‘Y'!V'V vL-I44'r4‘444l.!""'?I-c.44lv144o
’..‘!A\’\PQJeQPQxP1'PflQsRQIP4I... o4...4b441"viu.JOIOP'mueonOdO-Ordgme‘da
44...~9‘0o'PA'Da'-.\1P7,‘\c\V‘HQIVs 4.0.044644'64ud4nO.\o"o'¢'\‘4hd4ee
P'P"\o“"aefiti;nePof'Q‘Vfl'Veo. "2141‘ol.4‘4‘P’YP‘A4ctt'nroo4e.
P'P4"\!!DoonPi‘VntfiuR"e'PDPQpa «4.50.01.0"04Oc06544f9f:4oanesQQouw4Jc-e
ofl‘IRQQQQPP'3‘Vea§PP.‘I,1'!"'.O 4‘444.44h44l4ia'14b‘\b.vrt'4dchp
\l'AP'PD'JP"'D.Q‘I'IKPIQY‘P'VVQ .IKQCOIA4‘4JIl4-‘a\\b04st45441114.
\"'\‘¢P9'1u\'t'cYQVP'JPVI'HIIPK 4¢\Qi4rvaas44‘nd‘m444444414444ee
‘7"...'t“""‘\,l\“\fl.l.eauinn—R 5‘4445tICIQR44clhi4444'14444Iobit
"F7‘094)\,'eJ‘VVP’IQPVQODOQQJO.'r...’ 44iie4409afisOo4mualp-o.7990—44454.e
V'I‘QQI'VY‘QAVQIQ.e“'\l§eu.7.r'v bibsa”-075'!44ioa4P-!.Owio7...a.Irx4
F‘A'\"Plo,-"""V.'It'fi‘ntvelar 44-blce'Voe'99PeoD..o4aeper—-..
.A‘P‘lfl'a'Qoa5"!!!Qs'e'PPV-IVVV bui44c4‘pt..'4.t"v‘e.eoIVPI‘PQV..
"'\'3'Q¢.O.RVAD)-VIQ'ID-e'fl'n‘rfv 4s¢l41'0s'lOa'7—r‘QVQV-el\"ome...
"Y‘QRV'oOYI‘RCIOVVP.nPI'Ar\!I'v. e§4\.Oa\'et-n\‘roee‘npvap‘Yevrr.
"'\‘Rl“—"‘o‘497-I.'P"O‘Ya'Pp 00-49QmO4'P-“99s>\!II-e4o17v—11".
"Q.'\""‘QV§I’D‘QV;'VPIJCI'VY-a 4al.44vdfifr11vi‘Q‘Ve-D'O..vrwert-
‘Yfp'RVPIQAIe'Irnlnilhvo.DD'IQ-V .44o44~¢PpPQYRP4179uep.e.0woo»...
'QsYP‘VPRPe.'QVPQPVIVPp..aroe\.' 4—I—U.O.O.oJ.O‘.VVeOI.-.v1.V7...OOO.-..p'.
If)...PDQ.eeeugeeeNPCIP-P'Ievn'r- . er\....r.p-Ipenp'Qaqr.Os~...O~r.....

 

 

to 0.1 for

form T*

starting
and o=antivortex.

studies

O=vortex
(before quench).

quench,cfrom T=l.l).

Quenching
case.

K=0.1J

Fig. 3.5

quench).

(50 MCS/S after

(before
T=l.4).

T=0.l
o8

T=1.4>T

T=0.1 (50 MOS/S after quench, fr

T=1.1<T

 

 

(a)

45— Spin Direction

 

 

 

flea-

 

 

 

 

(b)
4-—- Spin Direction

 

 

Fig. 3.6 Two possible vortex-antivortex configurations.

56

57

-2.-

 

 

 

 

 

-4 .-
Ln(v)
-5 .-
'6 ,.
-7 1-
111(4)) -8 1 L A 1 A j
0 4 O 6 0.8 1 O I .2 1.4
0 [IT
1 ‘ Log of vortex-pair density vs II T to: 3600 spin
lattice. Straight line .3 least-squares fit to lowotempentute
\ dam.
5‘ 131.3
5 «1 °
4 1
3 d
O
2 d
0
1 <
O
. - . . . . . . e~=1/T
“5 10
Fig. 3.7 Log of the number u of VA pair vs l/T for the

32x32 spins system in K=0.1J case. The insertion is the
results of Tobochnik and Chester for isotropic XY model.

VI. HIGH TEMPERATURE SERIES EXPANSION

As an additional check, we have used the high temperature series
expansion (HTSE) to predict the critical temperature for our model in
the K=J limit. In this case, the Hamiltonian given by [3.3] reduces to
a simpler form

H 8 -(J+K)‘;)coseicosej, [3.16]

\o

and the calculations of HTSE becomes slightly easier.
For the partition function 2, according to the cumulant theorem, we

have

-Lan = J—ln i, dONexp[-fl(J+K) Z c050 c059 ]
N N -' a,» i j

I
-N—ln<exp[-3(J+K) (Epcosoicosojb

'3 1 <v2""‘>[12m [3.17]

W mu (2M)!

2x+

where

<V2>= <82) -<H>2

<V4>8 <84) -4<33><H> -3<Hz>2+12<32><fl>2-6<H>4

and asl/kBT. On a loose lattice such as square lattice, all the odd
terms vanish due to the correspondence between the ferromagnetic and
antiferromagnetic configurations. The values of <H2n> are then
calculated by counting all the graphs which contribute to the integral

and the weights associated with these graphs.(3'18)

The internal energy per spin U and the specific heat C can be

derived from [3.17] as

_ _ 3 - _ *3 232i _‘” 2n
U 7; 70.112). C 7; [3 3A -.§.anx ,

where x=fi(J+K).
Up to n=4, the specific heat is

c = 0.5::2 +0.8906x‘ +1.0970x6 +1.2822x8 [3.18]

If the phase transition is Ising-like, then we expect that the
specific heat has a leading term of the form(3°l7)
c . -Aln(l- xz/xé), (A>0), [3.19]
thus the coefficient an will have an asymptotic form
an::;:A/(nx%P)
or
ln(nan)-—»lna -2n lnxc= lnA +2n lnrg' ’ [3.20]

Using a2, a3, a given in [3.18], a linear regression gives

4

T; = 1.303.

Note this value is in good agreement with those given by MKRG and MCRG.

60

VII. SUMMARY

We hve studied the nature of the phase transition in a system
described by the Hamiltonian given in [3.3] with various theoretical
methods. The results obtained depict a consistent physical picture.
Our main results are as following.

The system possesses long-range order at low temperature. There is
only one phase transition in this model. Between the low-temperature
ferromagnetic phase and high-temperature paramagnetic phase, there is no
intermediate XY-like phase.

Besides the spin-wave excitation, the system has vortex-antivortex
pairs and domain-wall excitations. Near the critical point, the VA
pairs tend to get trapped in the walls suggesting aniattraction between
these two types of excitations.

The phase transition is Ising-like. This conclusion is based on the
facts a) the system iterates to an Ising-like system either in the
low-temperature ferromagnetic phase or in the high temperature
paramagnetic phase in the MKRG procedure; b)specific heat apparently
diverges at the critical temperature where the susceptibility also
apparently diverges; this behavior is not consistent with the
prescription of HT phase transition but is expected for an Ising-like
phase transition; c)the correlation length critical exponent 0
calculated by MCRG adopts a value close to the Ising value 1; d)Monte
Carlo quench study shows that the LRD is destroyed by the formation of
domains with opposite spin orientations as seen in the 2D Ising model.

MKRG and the energy-entropy argument suggest that if a

 

61

symmetry-breaking anisotropy, no matter how small it is, is present then
it dramatically changes the nature of the RT phase transition. In our
model, the phase transition becomes Ising-like-caused by the formation
of domain-walls whereas without the anisotropy it is a KT phase
transition caused by the unbinding of VA pairs. The ultimate reason for
this change is that in the presence of the anisotropy, the dominant term
in VA pair excitation energy is linearly proportional to the separation

instead of logarithmic in the thermodynamic limit.

62
CHAPTER 4
PHASE TRANSITIONS IN

DIATOMIC MOLECULAR MONOLAXER

I. INTRODUCTION

Phase transition from an orientationally disordered plastic phase to
an orientationally ordered ferroelastic phase in 3-dimensional (3D)
molecular solids occurs in many physical systems and is a well studied

phenomena.(4'l)

However, similar phenomena in 2D systems have only
recently been investigated in adsorbed molecular monolayers on graphite
substrate and the study of such transition is in infant stage. Typical

(4'2) and nitrogen(4‘3) on a graphite

systems are molecular oxygen
substrate. In fact the whole area of structural phase transition in
molecular overlayers is of great current interest.

One of the interesting characteristics of molecular solids is the
competition between direct and indirect (lattice mediated)
intermolecular interactions which lead to different types of
orientational ordering. Since the effect of this competition is

pronouced in 2D systems due to the enhanced fluctuations(4'4) on

(4

expects to see multiple stage transitions ’5) and fluctuation-driven

63

first order phase transitions(4'6)

in 2D systems.

In this chapter we report the results of the molecular dynamics (MD)
study of a 2D ferro-paraelastic phase _ transition by using
constant-pressure (force/length) ensemble. The system consists of 400
diatomic molecules whose mass centers and orientations are confined to a

2D plane denoted as the XY plane. The molecules interact through

atom-atom Lennard-Jones potentials,

vatom-atom= 4.[(a/x)12 -(a/R)6]. [4.1]

Adopting the parameters propriate to oxygen, this system closely
resembles the 8-phase of oxygen molecular monolayer absorbed on graphite
substrate. In the G-phase the axes of oxygen molecules are primarily
confined to the graphite substrate surface. The corrugation of the

substrate potential is known to be small(4‘7)

and has been neglected in
the present study.

In addition to the study of the nature (continuous or discontinuous;
one or two transitions) of 2D ferro-paraelastic phase transition we also
attempt to understand the nature of the topological excitations that
drive such a transition. Since the orientational modes of this system
can be represented approximately by a 2D anisotropic X!
Hamiltonian,("8) one expects, besides the librons and phonons,
vortex-antivortex (VA) pair, domain-walls and vacancy-like excitations
to be present and drive the phase transition depending on their relative
thermodynamic importance.

A promising candidate of 2D ferroelastic structure is the 5-phase of

02 on graphite. At low temperature, the molecules orient parallel to

64

each other and their centers of mass form an almost centered rectangular

lattice (see Fig. 4.1(a)) with b/a=2,45.(4°2)

which can be also treated
as a distorted triangular lattice. With increase in temperature it is
possible that the system undergoes a ferro-paraelastic phase transition.
In the paraelastic phase which can be denoted as an orientationally
disordered 2D plastic phase (we will use this name from now), the
centers of mass of molecules form an isotropic triangular lattice.
There are three possible ferroelastic structures which are related to
each other by rotating 2r/3 and the ground state is three-fold
degenerate.

If we describe the above transition by strain (the orientational
degrees of freedom being integrated over) then the strains

W = 314, - ). W = e .

l 2

283‘” KY

can be identified as the two components of the order parameter<4°9)

and the Ginzburg-Landau-Wilson Hamiltonian has the form

4 2 2 3 2 2 2
H6”, zrmlwil +—r2if111i+w<1,, 31, 12> +132.pr . [4.2]

This Hamiltonian is claimed to have a continuous transition belonging to

("9) Since for

the same univervality class as the 3-state Potts model.
the actual system consisting of molecules the above Hamiltonian is only
approximate and the configurations with antiferro ordering (herringbone
structure, (see Fig. 4.1(b)) ) have relatively low potential energy,
departures from the above predicted behavior are likely to happen, one
of the aims of this chapter is to explore this possiblity.

Furthermore, on qualitative grounds, one has reasons to conjecture

that the present system may show two transitions in going from the
ferroelastic phase to the plastic phase under suitable contiditions.
First the long-range orientational and ferroelastic order is lost by the
formation of domain walls between three possible~ ferroelastic domains,
in this phase the system has triangular structure on large length scale.
On a triangular lattice, the VA pair structure (in the effective X!
Hamiltonian description) is close to the herringbone (HB) structure,
(see Fig. 4.2) . Since the HB structure has an energy only 5.4% higher
than the ground state, or because of the presence of the phonon
excitations, the thermal average over all possible configurations of
domain walls generates infinitesimal relative rotations of nearest
neighbors in an average sense, it is possible that in this phase that
the vortices are still not free and form bound VA pairs and the
orientational correlation function only decays logarithmically. This
intermediate phase, which is similar to the XY-like phase present in

(4.10)

z(p) models when p>4, goes to the high temperature paraelastic

phase with exponentially decaying orientational correlation through a

hosterlitz-rhouless-like(4'11)

transition. Although the detection of
such a subtle phase directly in a MD simulation is difficult, we make
some effort to find if the ferro-paraelastic transition in this system
takes place in two stages.

The long controversial problem about the nature of 2D melting, i.e.,
whether the transition is first-order or second-order, is also studied
in this chapter.

In Sec. II we give the ground state structure and prove that the

long-range orientational order is present in our system. The molecular

dynamics simulation method we used is documented in Sec. III in great
detail. The existence of a plastic phase is suggested by the MD results
which are presented in Sec. IV, the ferro-paraelastic transition is
thought to be first-order. The driving mechanism for this transition is
dicussed. In the final section Sec. V, along with the results of the
simulations for the melting transition, which is believed to be
first-order, we make a few comments about the phase diagram of this 2D
diatomic molecular system.

In all the following calculations, the parameters we used are those
(4.?) .

reportedly appropriate for the oxygen molecule, i.e.,
. =54.34kB, 0 =3.055,
d 20.6045, m 815.995ma,
where ma is the atomic mass unit. When we talk about the

quadrupole-quadrupole interaction, for the computational convenience,
the quadrupole moment is constructed by putting ZZe charge on the
molecular center and -Ze on the two atomic positions. The value we used
is Ze=0.1ll esu, it corresponds to a quadrupole moment

c2=-o.39xlo"26 esu.cm2.

 

 

'[ - 67
4--------- } 1.208
1r,” \
I \
/ \
/ \
/ \
/ \
[LPN-b. ‘} \I
x 2 I
‘ l
\\ l/
\ I
.1. i---_§--_.? x
1— 3.332 --I
ammo mmnumn km
p=0.693

Fig. 4.1 (a) The gound state configuration. The_ t.of
energy U is c, the unit of density p is a . In
effective spin Hamiltonian, the spin orientation is the
twice of that of the molecule.

VA

|._ 3.975 --1

 

 

HEBRING BONE
E=-7.83

pm
Fig. 4.1 (b) The herringbone configuration which
corresponds to antiferromagnetic ordering in spin
language.

 

 

 

 

 

In spin representation, to form a VA pair in

Fig. 4.2 (a)
ferro-ordering state, at least one spin has to rotate

n/3.

 

 

 

 

 

Fig. 4.2 (b) In antiferro-ordering state, to form a VA
pair, spin 1 and spin 2 only need to deviate slightly

from the local minimum. The vortex, as topological
defect, is referred to ferromagnetic ordering.

n. _m_§ eaouwu sun:

(4'7) found that the

By using a pattern search program, Etters et al.
ground state of the low-density phase of oxygen monolayer on the
graphite substrate is a centered rectangular lattice (see Fig. 4.1(a))
with a=3.32A and b=8.07A. The oxygen molecules lie flat on the
substrate surface and are orientationally ordered along the b axis.
Because we believe our system is closely related to this 6-phase, we
start from this structure for searching the potential minimum. A nearly
zero value of the internal stress tensor is an additional criterion we
used for the location of the minimum. This procedure gives the same
structure as found by Btters et a1. with slightly different lattice
lengths, a=3.332A and b=8.054A. Calculated for a 19x19 lattice, the
potential energy is -8.28e with a surface tension of 0.0064e/02. The
internal stress tensor Pay and the surface tension p, (from now we will

call it "pressure”), are calculated through virial formula in its 2D

form,

P 3

up 7:? ( rmxfix" + 2 fix“), [4.3]

i i «,5: 13 13

1 . -’ .
and 9875(Pxx+Pyy)' where 0 is the area of the system, Fij ls the force
acting on molecule i by molecule j, 21351-3 j. If we add the
quadrupole-quadrupole interaction, the change of the lattice lengths is
less than 0.3%, and the change of the ground state energy is less than

2%. Taking the coverage of a V3x/3 structure as the coverage unit,

which corresponds 0.0636 molecules/A2 on graphite substrate, the ground

70

state of our system has a coverage of 1.172.

Using effective Hamiltonian, we found that the herring-bone
structure (see Fig. 4.1(b)) is energetically favorable for a triangular
lattice. We also searched whether there is a potential minimum for such
a configuration. The finding is that the energy has a minimum when
lattice spacing a is 4.0065, with a value of -7.38e, slightly higher
than the ground state energy.

Although true translational long-range order is absent in such a 20
system, it is possible for the system to possess orientational
long-range order. Using the harmonic expansion at sufficiently low
temperatures, it is easy to show that there is a gap between the ground

state energy and the libron excitation energy spectum. For example, if

the lattice is rigid, an effective spin Hamiltonian("8) for the
nearst-neighbor orientational interaction is
V(9i,Oj)=-0.388+0.625(coszoi+coszoj)-0.303(cos4oi+cos40j)
-0.584cos(291-20j) -0.409cos(291+29j), [4.4]

where Oi and oj are the angles between the molecular axes and the line
connecting the centers of these two molecules. In the previous two
chapters, we have already seen that such an anisotropic interaction
leads to orientational LRO. In the framework of the harmonic expansion,
if we take into account both the libron and phonon excitations, the

harmonic Hamiltonian can be written down as

Hharmg E €§¢E¢_§ + Z CE,§wfi¢§ + % e§¢§¢-§ ' [4'5]

where ibis the wave-vector of phonon modes and ‘3 is that of libron

71

modes. In the calculation of orientational correlation, the integration
. . 3.3

over wk in <8] > will leave a term exp(6ck'q¢q/4e2), so the effective

libron excitation spectrum can be thought as

e& = sq- E aficio'a, [4.6]

where “f is a coefficient, and ci,a.is expected to be linear to k in the
small x limit as in 30 case. The coupling with the translational
degrees of freedom can only reduce the gap with an amount proportional
to the square of coupling, except for some special case where the gap
may accidentally become zero. Thus if the orientational LED is present
for a fixed lattice, generally speaking, it will persist after we

release the mass centers of the molecules.

72

III. THE CONSTANT-PRESSURE MOLECULAR DYNAMICS

In comparison with Monte Carlo simulation, which is well documented,
in some aspects even seems standardized, and a few excellent review

books are already available,(4°12)

molecular dynamics has more room for
user's imagination and creativeness, some details scatter in the
literature or even remain in the research notes of experts. 0n the
another hand, I was witnessing a rapid growth of the applications of
computer simulation in this department in recent years, and further
growth in this tendency could be foreseen, as in every major
institution, when more and more computational facilities become

accessible to physists. For these reasons, I would like to give a ‘

detailed description of the MD method we have used.

A: Constant-Pressure Ensembles and Stress Tensor

for Multiatomic Molecular System

To allow our system to change its lattice structure as that occurs
in a ferro-paraelastic phase transition, we have applied a
constant-pressure molecular dynamics (MD) method first proposed by
Anderson(4'13) and later elaborated by Parrinello and Rahman(4’l4).
Take a monoatomic molecular system as an example, the main features of
this-method are as follows: first one defines a MD cell characterized by
vectors a and b, positions of the molecule are expressed in the basis
consisting of a and b instead of actual Cartesian coordinates. The

vectors a and b are allowed to change in the simulation, but at any

instant the positions of molecules are expressed as.

73

x ax bx (5(a))

(y) “ (aY by) 5(b) , [4.7]
where geax§+ayy and so on. Because of [4.7] if the MD cell changes
during the simulation, the positions of the molecule will be still well
defined. We will use the term "scaled coordinates” for s, the
components of s are in the interval (0,1). One then attaches a mass W
to the MD cell vectors. This is analogous to putting a piston on a gas
system, the motion of the piston is determined by the external force and
the pressure of the gas. To obtain the equations of motion, one then
needs a Lagrangian, and a possible one is

L = 1‘- 7. méfcé. -v + -1—WTrh'h -pa [4.8]

2.;, 1 l 2
where h is the 2x2 matrix appearing in [4.7], Gsh'h, 0 is area of the MD

cell, Dzlhl, and p is the external hydrostatic pressure.

The equations of motion generated from this Lagrangian are

3(8)
m'é.= h'1 z ( ‘3 -mG-lés. [4.9]
1 (W pf’?
13

This equation describes the motion of the i-th molecule with a
coordinate system whose basis vectors are a, b. For a and b themselves,

the equation of motion is

1

W1: 2 mew-pm' [4.10]

where

74

( 33ml}: if + z F'“.x".), [4.11]

p — ._1_
[JV [2 (3.5) i] i]

In eqn. [4.11] we formally write xi=hsi, ii=hsi. P v has the form of

u
the internal stress tensor calculated by. using virial theorem.
Eqn. [4.10] shows that the change of volume and shape of the MD cell is
driven by the imbalance between the external pressure and the internally
generated stress. The constant of motion is readily identified as the

enthalpy H given by
I mx.x.+v +p9 [4.12]
I

with a small correction ZkBT due to the 4 degrees of freedom of the MD
cell vectors.

The choice of MD cell mass W is theoretically arbitrary if only the
equilibrium properties of system are concerned because those properties
of the system are independent of the masses of its contitutients in the

classical statistical mechanics. As proved by Anderson,
average of a physical quantity along the trajectory produced by this
Lagrangian is equal to the isoenthalpic-isobaric ensemble average of
that quantity.

This novel simulation method has proved to be extremely powerful in
studying the phase transitions involving structural rearrangement as

beautifully shown by Rahman and coworkers.(4'l4)

al.(4°15)

Klein et

generalized it to the multiatomic molecular systems by adding

a term zwiIoi to the Lagrangian given in [4.8], where I is the moment of
t

inertia of molecules, oi is the angular velocity of i-th molecule. In '

75

this term, there is no coupling of the rotational dynamics with the MD
cell vectors unlike the term corresponding to the motion of the centers
of mass E siGsi, so [4.11] keeps the same form where only the motion of
the centers of mass contributes to the stress tensor Puv' It means that
no matter what kind of orientational order the system has, the volume
and the shape of the system will be decided directly by the motion of
the centers of mass. This may be reasonable for either a fluid or for a
solid which only undergoes an isotropic expansion or compression when
orientational order changes. For our system, where the disorder of the
molecular orientation is expected to play a crucial role in the
ferro-paraelastic phase transition, we realized the need to develop a
proper procedure to accommodate the contribution of the rotational
motion to the structural change.

For a rigid rotor, the rotation does not contribute to the
hydrostatic pressure, which therefore is determined totally by the
motion of centers of mass. But the rotation does make its mark on the
stress tensor. For our diatomic molecular system, if we write down the
equations of motion for every atom with the constraint of the rigid
rotor condition, then rewrite these equations in C.M. (R) and relative

(r) coordinates, for the latter, we have

0. A

- J- J -_D
mr- 2 (Pl 13‘2) +1yconstraint ' [4'13]

J .15
where F1 is the force acting on atom 1 and F2 on atom 2 of a particular
molecule. In polar coordinates, if we note that r=r=0 for a rigid rotor

and insert r=d, the above equation becomes

76

2

- ' e _L "’_ (U
”do 2 (F1 F2 )+Fconstraint’ [‘°l‘]
m9 —2 (F1 F2 ) . [4.15]

Eqn. [4.15] is simply the equation of motion for the orientational
variable 9, and eqn. [4.14] shows that in order to satisfy the rigid

rotor condition, the center of mass must provide a force with magnitude

- -mdéZ-L

_ (rL_¢r;
I"constraint 2(Fl F2 )' [4-16]

to atom 1 so the centrifugal force could be balanced. For atom 2 a
force with same magnitude but opposite direction must, be provided.
Taking account of the contribution of this force in the virial

expression for the internal stress tensor, we obtain

"0
fl

'0
I

(fl__(n
i1 P12)

E dCOSOiSln91(F

3 C.M. - (O) _ (0)
P P + § dcosais1n91(Fil F12)

cm. _ (O) _ 01
p . p a dsineisinoiw11 F32)

xY XY
= .n. m__¢u
Pyx ng + E dcoseiconoiw‘i1 F12) [4.17]

From eqn. [4.17] we can see that the hydrostatic pressure p which is
'%%(Pxx+Pyy) only includes the contribution of the centers of mass,

however Puv may have a significient part coming from the rotation. For

77

our system, at low temperatures, the molecular axes oscillate around
their ground state orientations, (see Fig. 4.1(a)), the term
-dcosoisinei(Fgg-Fgg) is always positive, so the contribution of the
rotation increases the xx-component meanwhile decreases the
yy-componenet of the stress tensor, this makes it easier for the system
to transform to a trianglar structure. To exclude this contribution in
simulations is equivalent to make the system artificially rigid, and
worst of all, if the coupling of the librons and the phonons reduces the
gap in the libron excitation spectum to a negative value making the
structure shown in Pig. 4.1(a) unstable under the perturbation of
librons, we might fail to detect this instability in the simulation
because the stress tensor was not properly calculated.

Suppose we connect the two atoms of a molecule with a very stiff
spring, so that the coupling of the orientational and translational (for
centers of mass) degrees of freedom with the vibrational degrees of

freedom is neglegible. In the limit the elastic constant of this spring

approaches infinity, we have a rigid rotor system described by the

Lagrangian
- 1 g, . 1 oz 1 0'. .-
L - 7Em51651+ 2 1;“. I91 v + 3mm h p0 [4.18]

with the internal stress given by [4.17]. The intramolecular vibration
contributes to the total energy a term ZNkBT which is omitted in our
calculations. Eqn. [4.15] and the equations of motion generated by
[4.18], i.e., [4.9], [4.10], [4.15], are used in our MD simulations

When the stress is calculated by using eqn. [4.17], the

 

78

. 1 .
antisymmetric component of the stress tensor o=-—(P -PYX) is zero for a

2 8?
system initially with zero total angular momentum as it should be. On
Cyber 750, through 104 time steps, the standard deviation of is of the

order 10-14

. This is probably the reason why some researchers who did
not include the rotatioal contribution to P“? found that the total
angular momentum was not conserved in their MD simulation in which the
Anderson-Parrinello—Rahman method was used to study the properties of
multiatomic molecular systems. For a system initially without net
translational momentum, the ”generalized momentum” g hsi should be zero

in the simulation, too. o and ; hs are monitored in our simulation in

i
order to check the accuracy of the algorithm used. It is important to
make sure that there is no net total momentum in the system. A non-zero
total momentum means that the center of mass of the system is moving,
this gives every molecule an additional velocity; it will increase the

pressure calculated by using virial formula and makes the system

artificially soft in the simulation.

8: Numerical Solution of the Equations of Motion

To integrate the equations of motion, we have used a

(4.16)

predictor-corrector algorithm whose accuracy is up to 0(At5) for

our case. If we denote the position variable of a molecule by qo, and

let

anAt)” cl" 0

 

where At is the time step used in the integration, then the predictor

values of the position variable and its derivatives up to order 5 are

79

given by

0 0 l 2 3 5 6
P 3 q *9 + q + q + q + q

pl= q1+2q2+3q3+4q5+ 5q6

p2= q2+3q3+6q5+10q6
p3= q3+4q5+10q6
p4= q5+ qu
p4= q6 [4.19]

These predictor values are input into the equations of motion to
calculate the accelarations, so we get the corrector value of the second

derivative c2. A=c2-p2 is then used to correct those predictor values,

cn= pn+ an [4.20]

where fn are 3/16, 251/360, 1, 11/18, 1/6, 1/60.

The advantage of this predictor-corrector method is that the
differences co-p0 and cl-p1 give estimates of the errors of the position
and the velocity. The procedure can be repeated if necessary, by using
the corrector as the next predictor. However, usually to predict and
correct once are satisfactory. When the external pressure is zero, so
that the enthalpy is equal to the internal energy, we found that in
104 time steps, the deviation of energy is less than 0.3%.

To initiate the integration, we need to know the initial positions
and their derivatives. To save computer time, it is desired that the
system is close to equilibrium even in the initial configuration.

Starting from the ground state configuration, we assign a small

positional deviation to every molecule. The deviation can be a product

of a fixed value and a random number with approximate Gaussian
distribution. The fixed value is time interval At multipled by the mean
square root velocity at required temperature T0 according to the
Maxwell-Boltzmann distribution. The time step used in this “initiation
stage" is one tenth of that used in the simulations. The approximate

(4.17) The

Gaussian random number is generated by a simple procedure.
central limit theorem says that the sums of n random numbers, no matter
what kind distribution those random numbers have, when n approaches
infinity, will reach a normal distribution. Using the available
function subroutine RANP( ) on Cyber 750, which generates a random
number in interval (0.1) with uniform distribution, we take a sum of 12
such random numbers. Since the expectation and variance of the uniform
distribution on interval (0,1) are 1/2 and 1/12, the sums will have
expectation and variance 6 and 1. By subtracting 6 from the sum, we
have an approximate standard normal distribution, almost
indistinguishable from a true Gaussian distribution by eye, with the
desired feature: no number has a magnitude bigger than 3, no molecule
will have a initial velocity larger than three times the mean square
root velocity. The molecules now move to the new position q1 form the
ground-state position qo, the velocity is given by
v = ( ql-q0)/At (q=x,y,9). [4.21]

To make the total momentum zero, we shift all the velocities, i.e.,

ii

1
v.= v ——— 2 v
llN;

where N is the number of the molecules. Since the velocities have an

81

almost Gaussian distribution, 2 vi is close to zero, and therefore the
1.

shift is very small. Then, the velocity vx is multiplied by a factor

2 Wviyxi/ Z vixyi to make the total angular momentum zero, Since

2 mdzé is already adjusted to zero. Finally, the adjusted velocities

i
are used to calculate the temperature T of the system by the
equipartition theorem
l
__ __ z '2
:NkBTs 2 §m(vfx+viy+d219. ), [4-22]

All the velocities then are multiplied by a factor VT07T to force the

system to have the required temperature T , which is not too far from

0
zero.
A simple integration procedure, namely
qn+1= 2qn'qn-l +qn
can be used with the small time step At to obtain the first few

configurations. Then these positions are transformed to position and

derivatives at time step n by the equations

 

 

 

 

 

go o o 120 o o o qn+2

q1 -6 so 40 -120 30 —4 qn+1

q2 1 -s so -150 80 -s o qn

q3 ““12; s 5 -50 70 -3s 5 qn_1

q‘ 5 -2o 30 -20 s o qn_2

qs 1 -5 10 -lo 5 -l qn_3 [4.23]

Note that the superscript indicates the order of the derivative, the

subscript labels the time step. At this stage, the calculation is

82

performed with a fixed MD cell corresponding to the ground state
structure. After obtaining the high-order derivatives, we change the
value of At in q“ from the small initialvalue to that used in the
simulations. This finishes the "initiation stage" and the equations of

motion are now ready to be integrated.

C: Aging and Generation of the Equilibrium Ensembles

The next stage of simulation is the "aging” stage as called by
Rahman. The volume and the shape of the MD cell are now allowed to
change by introducing the mass W attached to the MD cell vectors. W
will affect the relaxation time of fluctuations of the MD cell, and in
our simulations we choose W- min. This choice is found to be reasonable
for our system because in simulations, the vibration period of the MD
cell is smaller than the length of a typical MD run (2000 time steps),
but much larger than the correlation time of the dynamics of molecules.

In the aging stage, after every m time steps, where m is a
predetermined number, the instantaneous temperature of the system is
calculated by using [4.22], then all the derivatives are multipled by a
factor of (5375 to force the temperature to be To. Because the system
is initially in a state not far from the ground state, the temperature
may rise rapidly at the beginning while a considerable amount of
potential energy transforms to kinetic energy. The procedure of
adjusting all the derivatives instead of the velocities only, shortens
the relaxation time of this tumultuous period, and helps the system to
reach equilibrium much faster. After running the system for a few

thousand time steps, we switch off the temperature adjustment, and run

83

the system for several segments (each segment consisting of 500-1000
time steps) to see whether the averages of energy, temperature, pressure
are stable and the kinetic energy is equipartitioned. If the result is
unsatisfactory, the aging process is repeated. At a temperature not too
close to the transition temperature, aging the system once appears
sufficient. However, near the transition temperature, for example, in
the solid-liquid coexistence region, 104 time steps may be needed for
the system to melt to liquid from a solid state. After the aging
process, the system is used to calculate the equilibrium or dynamic
properties under study.

When the system is heated or cooled to a new temperature, the aging
process is carried out again at this new temperature. Usually, except
for quick quench studies, we give SOD-1000 steps for the purpose of
raising the temperature. These steps are divided into 10-20 segments,
the temperature difference is uniformly distributed among these
segments, so the heating (cooling) process is smooth and gentle. The
cost of a rapid heating (cooling) is always an aging process intolerably
long, and even spurious results in some cases. Close to the
transition,the temperature step we used is as small as 0.02.

The final stage which we call as the "equilibrium stage", where the
fluctuation of temperature is small. After m time steps, we adjust the

instantaneous temperature to the desired value T and record the

0
information needed on tape for later calculations of the
isobaric-isothermal ensemble averages. Usually, the average is
calculated over 2000-3000 time steps, however, close to the transition

temperature, 1.0-1.5x104 time steps are used.

84

D: Simulations of the Lennard-Jones System
Following the convention in MD simulations of the Lennard-Jones

(4°14'4'18) we use the reduced dimensionless units for

systems,
computational convenience. It also has the advantage that the
simulation results are independent of certain parameters, for example,
the energy parameter s and mass m. The lengths are expressed in units
a, the energies in units e, the temperature in units e/kB and pressure
in units c/az, m is the mass unit, these choices lead to a time unit
an/e, which is 1.8148xlO—lzsec.. The high accuracy of the algorithm we
have used and the relatively large cut-off distance 50 (a popular value
is 3.20) for the interaction potential enable us to use a relatively
large time step, At=0.01102, corresponding to 0.02 picosecond, even when
the fast librational motion is present.

The Lennard-Jones potential V(r) is cut off at r =50 and then

0
modified to be of the form

, V(r)-V(r )+f(r )r-f(r )r , r<r
V(r)={ o o o o o

0, r>r [4.24]

0

so that the potential and the force are both continuous at the cut-off,

where

avm
ar

 

f(r)= -

Since f(ro) is as small as -0.0003l, the term linear in r in [4.24] does

not affect the calculation significiently. In fact, this effect is

expected to be smaller than the statistic error, but the continuities of
potential and its derivative are helpful to the accuracy of the
integration. The large cut-off also gives more information about the
radial distribution function of the mass centers (RDch), because the
RDFcm is calculated along with the potential and force for the reason of
computational convenience when pair-wise potential is used.

The major amount of computer time is spent on calculating the
forces. To reduce the computation time, we divided the system into 16
blocks, the length of the block in one of the two directions (a and b)
is five lattice spacing, which is much larger than the cut-off distance.
The molecules in one block can only interact with the molecules in the
same block or the neighboring blocks. The molecules are labelled from 1
to 400 and are distributed in the blocks according to their positions.
In the central memory .of the computer, every block has 40 addresses,
this is sufficient to accommodate the labels of the molecules in a
block, because the average number of molecules in one block is 25, with
a possible variance of 5. In the central memory, these labels are the
contents of the addresses symbolized by the block name. In
calculations, the particular pair of molecules which interact with each
other can be easily identified without searching all the molecules and
comparing the separation with cut-off distance 5 for every pair of
them. Thus the number of loops in calculation of the potential and
force is proportional to N rather than to N2.

When a molecule moves out of the MD cell, it should be replaced by
its image according to the periodic boundary condition. However, to

calculate the correlation functions of position, we have to know the

positions of every molecule even some of them may walk out of the MD
cell. Rather than to add an array to count the number of
transgressions, we exploit a hardware feature of the Cyber 750. A word
in Cyber 750 has 60 bits, for a floating point real number, the first is
a sign bit, the subsequent eleven are used to store the exponent and the
remaining 48 store the coefficient. If we shift the the position
variables (in the scaled coordinate) to be in the interval (24, 25), and
the angular variables can be scaled in the units 2w, then shifted to the
interval (24, 25). too, thus the first 5 bits of the 48 bits storing the
coefficient will contain the integer part of the position variable.
Even for a very long simulation run, it seems unlikely that the integer
part can become less than 16 or larger than 32 resulting in a change of
the number of the bits storing it, because a length of 8 in our scaled
coordinate system corresponds to 160 lattice spacing! In fact, we have
never found it happen in our simulations. Knowing exactly where the
decimal parts are stored, it is trivial to fetch them only into the
calculation of potential and force. Now the decimal part takes care of
the periodic boundary condition, the integer part keeps the information
whether the molecule flees out, altogether, they give the true position
of the molecule. These data can be then used directly to calculate
position correlation function without further manipulations. The only
compromise we have to make is that the 14 significient digits Cyber 750
provides are now reduced to a number a little worse, but still better

than enough, 12.

87

IV. THE PERRO-PARAELASTIC TRANSITION

Under zero external pressure, beginning at T=0.12, the temperature
of the system is gradually raised with a temperature step 0.06. We
found two phase transitions, the first one at T1-0.38 (20.6K) where the
system undergoes a ferro-paraelastic transition, the orientational LRO
is lost and the lattice becomes triangular; the second transition

occured at T =0.70 (38.0K) where the plastic phase melts to a 2d liquid.

2

Various physical quantities such as internal energy, density, radial
distribution function, order parameters are calculated in the MD
simulations.. For the ferro-paraelastic transition, the results are
consistent with a first-order behaviour, and do not show any clue of a
two-stage transition. When the system is heated from T=0.36 to. 0.38,
all these quantities changes drastically in 8000 time steps, this
simulation length corresponds to a real time 160 picoseconds. After
that the system has a triangular lattice and the orientational order

disappears. The system then run 1.2x104 time steps over which thermal

averages are calculated.

A: Lattice Anisotropy Parameter and Orientational Order Parameter

We use 7=(b//3a-l) as a measure of the distortion of the lattice
from a perfect triangular structure. In the ground state, a=3.33ZA and
b=8.054A, 7 has the value 0.396. Just before the transition, the MD
cell gives a=3.462, b=8.055 and 7 is as high as 0.343. The slow
reduction of 7 as T increases (see Fig. 4.3) proves that the ground

state structure is stable under perturbations of librons and phonons.

After the transition, 3 increases suddenly to 4.006 and b is 7.066,
these values correspond to a triangular lattice. An interesting thing
is that the lattice spacing 4.006 is very close to that of the HB
structure in Fig. 4.1(b). The curve of 7 vs. T is plotted in Fig. 4.3.

For the orientational order, we take

1

722: —N. E: (105(291-29) [4.25]
as the primary order parameter, where 5=ޤ I 01, and
I
n = 1‘- }: cos(69 455) [4 26]
6 Al i i . °

as an auxilary order parameter. If the orientational LRO is destroyed
by the formation of the three equivalent domains and the
ferro-paraelastic transition takes place in two stages, our test
simulation studies show that "6 has a sharp but relatively smaller drop -
at a critical temperature as n2 does when the domain walls appears, it
will decrease quite slowly until a higher critical temperature where it
drops sharply again as the structure of the system becomes triangular.
In this intermediate phase, n2 keeps a value close to zero.

For the present system, both n2 and "6 have large drops at T=0.38,
it contradicts the hypothesis of a two-stage transition. In Fig. 4.3
the orientational order parameters appear discontinuous at the
transition temperature 0.38, and n2 does not show a slow decaying tail.
It is a well-known phenomenon that such a tail is always present after
the transition as a finite size effect in computer simulation if the
transition is second order.(4'12) These facts suggest that the

transition is first order. The tail of "6 is due to the local short

range ferro-ordering along three equivalent directions. After
transition, in the system there are small patches where this ordering

persists.

B: Radial Distribution ruction

Because the ground state structure can be treated as distorted
triangular lattice, the radial distribution function of the mass centers
of molecules (RDPcm) can be viewed .as that every peak in a set
characterizing the triangular lattice is split to a doublet. In
particular, the peak coming from the 6 nearest-neighbors in a triangular
lattice will split into two separate peaks, one coming from two
neighbors with separation a, the other from the four neighbors with a
distance (a2+b2/4)1/2. In the case where the formation of domain walls
destroys the orientational LRO, the radial distribution function of the
mass centers of molecules (RDch) will have two sets of coexisting
peaks: one reflecting the distorted triangular lattice structure (coming
from the domains) and the other characterizing the triangular structure
(coming from the wall region). Especially, between the two peaks
characterizing the 6 closest neighbors in the distorted triangular
lattice, emerges a new peak signaling the appearce of the triangular
structure in the wall region. This is what happens in some of our test
runs.

In the first 8000 time steps when the system is heated to T=0.38
from 0.36, the splittings between the two close peaks which correspond

to a single peak of triangular lattice becomes narrower as time

increases, finally the splittings disappear and the RDFcm shows a set of

peaks appropriate to a triangular lattice. Fig. 4.4(a) and (b) give the
RDFcm at T=0.36 and T=0.38, respectively

The explicit form of the radial distribution function g(r) in

Fig. 4.4 is
<n >= «Arzg(r) [4 27]
r2,r3+Ar2 .
where nr r +Ar is the average number of molecules in the ring whose
I

area is zArz. Because the horizonal axis is plotted in r2, g(r) will

approach the value of the density as r2 goes to infinity.

C: Energy and Density

In Fig. 4.5 and Fig. 4.6, the curves of energy E and density p vs. T
are given. Both show a jump at T=0.38. The change of the entropy can
be calculated by using the formula AS=AE/Tc and it is 0.88 (in unit k3).
Since for a free rotor the total entropy is ln(2r)-l.83, the
considerable change of the entropy suggests that the rotors are nearly
free. This fact, along with the discontinuities in the energy and the

density, prefers the conclution that the transition is first order.

D: Orientational Diffusion Coefficient

In the ferroelastic phase, due to the symmetry of the diatomic
molecule, there are two degenerate orientational states separated by x
for every molecule. As the temperature increases, so does the
probability for the molecules to overcome the energy barrier and shuttle
between these two states. Before the rotors become free, one expects

them to behave like hindered rotors and show a small orientational

91

self-diffusion coefficient at large time scale at certain temperatures.
If the rotors are free, then even at a relatively shorter time scale the
system will have a considerable orientational self-diffusion
coefficient. '

The orientational self-diffusion coefficient D0 is expressed in the

well-know Einstein's relation (in its 1D form)

D9=lim441— <[9(t)-9(0)]2>. [4.28]
t+~* t

In MD simulation, the average of orientational autocorrelation function

is calculated as

f(t)=<[9(t)-9(O)]2>8”.Z.‘—1- Z ‘1' [9 (t+kT)-9 0(7)]2 [4 29]
‘fl,flu[i I“ i i ' °

where r is a constant decided by the correlation time of the molecule
dynamics, mk is limited by the number of configurations available. Then

we use linear regression to derive D from f(t) with different time

9
scales.

When T is 0.30, at a time scale as large as 1600 time steps (32 ps),
there is no evidence of the occurrence of orientational diffusion.
Heated to T=0.36, we found D9=4.910.4x109rad2/sec for t>1000 time steps,
at T=0.38, D9 jumps to 1.49:0.08x1011rad2/sec after the same number of

time steps (see Fig. 4.7). Fig. 4.8 gives lnD vs. 1/T, note the slope

8
changes at the solid-liquid transition region about T=0.70.

E: Strain Fluctuation and Elastic Constant

The gauge metrix G appearing in eqn. [4.8] reflects the change of

the volume and shape of the MD cell, it keeps the information about the

 

92

strain. As Parrinello and Rahman have identified, the relation between

G and the strain tensor can be expressed as

f(t)=—[hb-lG(t)h31 1] ' [4.30]

where ho characterizes a reference MD cell, Rahman has recommended to

“59 the average of h(t) over the simulation as h0.(4'14)

The compliances are related to the fluctuations of the strain by the

fluctuation-dissipation relation as(4‘19)

I)

r. . =_<‘. .8 > [4e31]
iJkl ‘6T .13 kl

and the elastic constant C can be obtained from the inverse of

ijkl
rijkl' However, the 4x4 matrix of rijkl is Singular due to the
permutation symmetry for the Cartesian index, by using the Voigt
convention(4'2°)

r ={r1jkl' 1(8). (052)
pa 3
zrijkl, k 1 (0)2) [4.32]

one can get a 3x3 regular matrix and its inverse can be calculated. To

obtain symmetric C one then inverts the transformation of [4.32].

ijkl

To derive accurate values of Cijkl

fluctuations of strain near the transition. However, the general

is difficult due to the large

features can be seen by comparing the values at T=0.36 and 0.38. We

found that (in unit e/az, in terms of Lennard-Jones potential

parameters)

 

93

for T=0.36, C =l.45, C

11 =0.97, C

=2.lO, C 481.33;

22 12 4

for T30.38, C =1.05, C =1.07, C =1.02, C4‘80.11.

11 22 12
Clearly, before the transition the system is' highly anisotropic but
after the transition it is isotropic, Cllzczz, but the ratio c44/Cll and
(Cu-Clzflc11 drop sharply at the critical point. This fact can be
visibly recognized from Fig. 4.9(b) where we have plotted the bonds
connecting neighboring molecules, the amplitudes of phonons are much

larger than that in Pig. 4.9(a).

P: Discussions
The prediction of the existence of a plastic phase from our MD
simulations is consistent with the Monte Carlo (MC) simulation results

of Etters et al.(4'7)

for small clusters. In the simulation studies of
clusters with different numbers of molecules, they found that there was
an orientationally disordered phase when the number of molecules in the
cluster ranged from 3 to 12 (the largest number in their simulations).
Our MD results and their MC results are complementary to each other:
they mimicked the behavior of the cluster using a more complicated
02-02 potential which included both a Lennard-Jones part as ours and a
small quadrupole-quadrupole part, the absorbate-substrate potential was
also added in their simulations; we simulated the monolayer including
only the dominant part of the 02-02 potential. For the clusters, the
the orientationally disordered behavior begins at about S-7K for the

molecules having 2 neighbors, and this temperature rises to ll-l7K for

molecules connecting 3 or 4 neghbors. For the monolayer, our simulation

94

shows that the transition takes place between 19.6 and 20.6K. The
question whether this plastic phase exists in real oxygen monolayer
remains yet to be observed experimentally. An interesting feature of
bulk 1-02 worth mentioning here is that an orientational transition is
believed to be present in this material prior to the melting temperature
54.8K.(4'21)

The discontinuities found in energy, density, orientational order
and lattice distortion parameter, orientational diffusion coefficient
suggest that the transition is first order. considering the evidence
that in 2D systems fluctuations may drive the orientational transition
first order for a fixed lattice (a system believed to represent this
situation is the N2 monolayer physisorbed on graphite substrate(‘°6)),

coupling to phonon excitations is more likely to make the transition
first order due to enhanced fluctuations.

To investigate the physical reasons for the observed first order
transition we have counted the vortex density on triangular plaquetes as
a function of T. Before the transition i.e. at T=0.36, there is only
one VA pair but after the transition i.e. at T=0.38 the number suddenly
jumps to 83. (see Fig. 4.9) This phenomenon can be understood as
follows.

In the low T ferroelastic configuration on a triangle, if one
molecule is fixed, then to form a vortex, one of the other two have to
rotate at least n/6; but in a herringbone structure, if the molecule
which has an orientation different from that of the other two is fixed,

then the other two only need to rotate a small angle but in opposite

directions in order to form a vortex (see Fig. 4.2). We have calculated

the potential wells in above two situations, assuming the neighboring
molecules are fixed (see Fig. 4.10). The potential well for the
ferro-ordering on the centered rectangular lattice is much narrower and
deeper when the lattice spacings are close to the values appropriate for
the ground state, but at T=0.36, it has nearly the same width as that
for the herringbone ordering on the triangular lattice, even a rotation
of «l6 does not cost too much energy. Therefore at low T the VA pairs
are strongly suppressed, without enough sources of strings, it is

difficult for the formation of domain walls.(4°22)

As T rises, the
orientational potential well broadens and becomes shallower and this
change is not only due to the orientational fluctuation, but also
drastically depends on the increase in lattice spacing caused by
anharmonic phonons. At a critical temperature TC, the core energy of
the vn pair is sufficiently small so that the VA pairs appear in a large
number and cause a first order transition. The .coupling between
orientational and translational degrees of freedom effectively reduces
the core-energy of VA pair and drives the orientational order-disorder
transition first order.

The low core-energy vortices tend to form clusters as can be seen
from Fig. 4.9(b) and as observed by Saito in his dislocation vector

system.(4’23)

Thus in many areas in the system the local density of VA
pair is very‘ high. For ferro-ordering, the vortex, as topological
defects, has a configuration which resembles closely the herringbone
structure, and on a triangular lattice, this structure is the ground

state for antiferromagnetic interaction. By coupling to the

translational degrees of freedom, the nature of the orientational

interaction can be changed as the lattice spacing varies. The
intermolecular interaction prefer the ferro-ordering at relatively small
separations but favors antiferro-ordering at. large separations, the
crossover is at about 4.0A. We have already given an effective spin
Hamiltonian for a separation 3.332A in eqn. [4.3] and it is
ferromagnetic in the spin language. When the separation is 4.03, the
effective Hamiltonian becomes

V(ei,9j)= -l.826+0.722(co5291+co529j)+0.07l(cos4oi+cos49j)

+0.3S3cos(201-20j)+0.569cos(201+26j) [4.33]
and it is antiferromagnetic. Thus in the region where the local density
of VA pair is high, to reduce the excitation energy it is favorable to
change the local structure to be triangular with a lattice spacing of
about 43. The net result of this lattice-mediated competition between
ferro- and antiferro-ordering is that before a possible 3-state
continuous Potts transition can take place, the ferroelastic phase
undergoes an orientational order-disorder transition and becomes
paraelastic simultaneously. The ferro-paraelastic transition becomes a
synonym- of the orientational order-disorder transition in this
particular case in the sense that they are completed at the same
critical temperature.

A configuration of a quench study is given in Fig. 4.11. The system
is quenched from T=0.38 to 0.01 in 20 time steps by taking the kinatic
energy out. There are only tiny domains with three equivalent
distortion directions; but a great number of vortices is persistent in
the system, and we find a hexagon with local herringbone structure. The

picture is quite different from that appearing in Fig. 3.5.

97

 

 

 

 

 

 

 

 

k; - 7 . 0.4
1.0 \
\ \
\.\. 772
B \
0.8 - \ o - 0.3
\
‘\
\ O
UJIr “\. ‘7
b. a
7;. \ 0.2
\O
OJI- \
\
- 00'
0.2 - \
D.°\
~o ‘ ~o ~ _
—' O
. on I l i
ne'o o.‘ 0e2 0.3 0.4 0.5 0.6
T‘
Fig. 4.3 Lattice anisotropy parameter 7 and orientational

parameters n2 and "6 vs. temperature. All parameters
sharply drop at T=0.36.

 

$1.00 an sine 11.“ up.» 1pm 1:00 1?.“ 29.00 :3.“

.4.

Reduced Teeperattre ie 0.3
n Steps A»: From iaeei To issue

amulet 13L""§9Jo‘°“ qugf'rtor‘m

 

 

 

 

 

 

 

I
8o”

8
0o”

 

 

. 4m {.se e‘.ee 1'." b. h . . .
[NR-81.“ flies 17-60 flee rise

Fig. 4.4 (a) The radial distribution function of the mass
centers of molecules at T=0.36, averaged over 3000 time
steps (60 ps). The positions of the peaks indicate a
centered rectangular lattice.

gum else s‘.ee '1.“ up.» 1:50 15.00 13.“ q.” 2;.“

 

 

 

 

 

8
.q ’2
Reduced Tenet-attire 1a 0.3
8
3‘1 Steps We From 1M1 To 1” :3
H
E.»
24% [‘8‘
z W
53
'53.. 3
E?' r;
1
Si
.1 P“
c
E:
93.1 '8'
B
61 ,8;
a J z
‘Imo '.

 

sise sloe 1'.“ than t . . .
R 1.50: than 1')“ shoe ace

Fig. 4.4 (b) The radial distribution function of the mass
centers of molecules at T=0.38, averaged over 3000 time
steps (60 ps). The positions of the peaks indicate a
triangular lattice.

100

 

 

 

 

-3.c . _ 1-
.
4.04
5.134
-am
4.0 ‘
-em
Eu
Fig. 4.5 Total energy E vs. temperature. Energy and

temperature are in unit e and «s/kB respectively. At
T=0.36 there are two data points, one is the data of
heating run, the other is the data of a cooling run
(cooled from T=0.38). The line are guides to the eye.

-

 

0.70

 

101

 

 

0.804
’—
0.” f V T I v v t
0.0 u 0.: 9.3 0.4 as cs 0.1 of. is 1.0
Fig. 4.6 Density p vs. _ temperature. Density and

temperature are in unit a and e/ respectively. At
T=0.36 there are two data points, one is the data of
heating run, the other is the data of a cooling run
(cooled from T=0.38). The line are guides to the eye.

(RIENTRTIM. MWTIW WIN
8.42

 

102

0.x -

0.x -

8.3 ‘
Ea d

0.3

0.3 "

 

 

 

148% 1689 1889 229a ammo

Time Step

 

11‘

 

 

 

14e3 1608 1893 amen 2299

'Time Step

Fig. 4.7 Orientational autocorrelation function from 1400
to 5300 time steps. The vertical axes are in unit of
rad . The temperature for the graph at top is 0.36 and
it is 0.38 for the graph at bottom.

103

 

Ills

lllll

J_l

 

 

 

1.0 2.0 _ 3.0

Fig. 4.8 Orientational diffusion coefficient D vs. the
reciprocal of temperature. Note D is discogtinuous at
solid-liquid transition temperature 8.7.

104

¢'+III‘77hqhdihq¢7I7III<
I'II-IIII‘IIqt¢¢qq¢rf¢\g¢7.h
II7IIII7HA¢¢7V¢¢¢I7777<
77777777 I¢q¢¢¢7IIIh7§7
777II<77wdntdqq¢7777¢77
IIIYIIIII77I87IIIIIIII
7577777777I777777I777
77777777717777.77777777777
IITIIIII7¢7V+III‘LII‘I'IIEI
(«T-IIIITI'I’I‘77"{£'IIIIIII
I; I I W I T I I I I I I I‘ ] I T 7 7 Y I 7‘ 7
my. :7. r7 7 7 I I ‘7 7I <7- 7 I I I I *7. 7 7 I «1-
7 7 7 I I I 7 C \ 7: I 7- 7 7" I 71 7' 7- 7 7' 77
41777-777”7777777I4I-777I'
I? ’3: 7 q- I" I I I I 7' 17> I I 7*- I‘ I' I I 3"" I I (I
‘ .' '7, 7 7 I 7 «7 7 rs; 7: I 77
r] 7 : 7 7' 7 I’ I I- 7: q- I <7 <7 I" I I' 7' I I I I'
I“ 1") I I I 7‘ V' ‘7' ‘T‘ 7 7 I =7 7

 

I

«4

e.
’v-r
M
.57:-
x’
.c;
”r
_.

41)
~9
H)
.4' ~
- ‘ ', ,
_.."+
49‘

 

¢L.§ 7 z a i k y A. T ;. f. f f I q 7 c
¢ & I‘ I’ 7 I Q Q Q q ¢ ¢ I I ¢ V I N. 7 '7 ¢ h

L

I...

Fig. 4.9 (a) A snapshot of the configuration of molecules
at T=0.36. The time step is 15000. In the picture,
circle denotes positive vortex and triangle denotes
negtive vortex, there is only one bound VA pair. The
frame indicates the MD cell, the molecules outside of the
frame are plotted according to the periodic boundary
condition.

¢Ia¢¢a”I*“””””¢I

 

. AA A' I “6.,

. 45mm!” vmvmmm
"1‘1 , ‘ AV” WAVAVVVfif

.7 ANAVAV‘AV .m an“:

     

Fig. 4.9 (b) A snapshot of the configuration of molecules

at T=0.38. The time step is 12000. After
ferro-paraelastic transition, the number of VA pair
suddenly jumps. In this picture, there are 83 such

pairs.

105

a Z antp

  
 

 

 

 

 

15¢ "‘
1

106 -
2

56+
// \‘3

/
/
/

cc.

0' 36' sh‘ 90'

Fig. 4.10

Orientational potential well for molecule.

Line

1 corresponds to the ground state, line 2 to the centered

rectagular lattice at T=0.36,
(lattice spacing is

triagular lattice
antiferro-ordering.

dashed line 3 to the
4.006) with

106

107

 

 

 

 

 

' 39 ’ o A? ' . I. "
‘0' (b 1‘ ‘1‘040 6 5? w "& 2’ Q ¢ ‘9‘ E ’2 .e' (P ? ‘St

RR¢¢3552Q¢9r¢§9~®559 gfiu
Q £A¢£2>¢ é; ¢ ¢ ¢ i! z'.z'13 2r ¢ ¢ 3‘ ¢ s.-g_‘3 Q Fléfl ;>'
‘ $5 ' ’1‘ ,3

-, . . . ‘6— 3 5 ¢ 3 3,3 2’ gfl¢~&~\&.&.grx&o—QJIJ‘

Fig. 4.11 A configuration after the system was quenced form
T=0.38 to 0.01. Although 1000 time steps are passed, the
density of VA pair is high, and the domains are small. A
hexagon with local herringbone structure is marked by
dashed lines.

108

VI. SUMMAR! AND COMMENTS

We have also used the constant-pressure molecular dynamics method to
explore the melting transition. We have reached the liquid phase at
T=0.70, where the RDFcm has less features, (see Fig. 4.12). In
Fig. 4.13, we can find many dislocation pairs in the configuration of
molecules. Fig. 4.6 and 4.7 show that the energy E and the density p
both have a discontinuity at the melting temperature, these results
indicate that the melting transition is first order as is seen in the
computer simulations for 2D monoatomic molecular systems with

Lennard-Jones potential.(4’24)

The time scale in our simulation,
however, is not sufficiently long to give the dislocation-mediated
melting theory an exclusive test.

According to the results of our constant-pressure molecular
simulations, the 20 diatomic molecular system studied first udergoes a
ferro-paraelastic transition from an orientationally ordered
ferroelastic phase to an orientationally disordered plastic phase at
temperature Tl, then at a higher temperature T2 the plastic phase melts
to liquid. Both transitions are first order.

The driving mechanism for the ferro-paraelastic phase can be
understood in the light of the role played by vortices. The sudden
emergence of great number of vortices destroys the orientational
long-range order. The physical reason why the transition is first order
and why it does not take place in two stages can be found in the strong

coupling between the orientational and translational degrees of freedom.

This coupling reduces the chemical potential needed to generate VA pair,

109

so the orientational order-disorder transition occurs at a temperature
much lower than that for a fixed lattice but with the same kind of
molecular interaction(4'25) and the transition becomes first order.
This coupling also rules out the possiblity of a two-stage transition by
introducting the lattice-mediated competition between different
orderings.

The reason for this transition to be first order is somewhat similar
to that proposed by Van Himbergen(4'26) for an isotropic XY model with

an interaction of the form
- 2 q2
V(ai-9j)-V(e)=2J[1-(cos 9/2) ], J>O [4.34]

for q2>lo in which case the transition from RT phase to the high T phase
is first order. For q2<10 , the potential is rather soft and the RT
transition is continuous. It will be very interesting to investigate
the effect of the symmetry-breaking field terms of the form hpcos(pei)
on this Hamiltonian.(4'27) Whether a large q value will change the
continuous Ising-like transition. to become first—order when p=2, or
whether at certain q the intermediate XY-like phase will disappear when
p>4? These crossover behaviors between different types of phase
transitions are of current theoretical interest.(4'28)

For diatomic molecules interacting with each other via atom-atom
potentials, the interaction is highly anisotropic when the average
separation between the molecules is small, but the anisotropic part of
the interaction reduces rapidly as the separation increases. The

relative strength of the anisotropic interaction can be measured by the

ratio between the magnitude of of anisotropic terms and that of the

110

constant term in the effective spin Hamiltonian. This ratio for
eqn. [4.33] is much smaller than that for eqn. [4.3]. Under zero
external pressure, the system can expand freely as the temperature
rises, and a tiny increase in the lattice spacing is sufficient to bring
the interaction to become close to isotropic, but the translational
potential well is still quite deep in comparison to the kinetic energy,
so the topological defects mainly appear in the orientational space, and
it leads to a stable isotropic plastic phase between anisotropic solid
and isotropic liquid phase. However, if the external preSsure (or an
uniaxial stress along a direction) is sufficiently high, the system is
considerably compressed from its zero-pressure state, the interaction
will keep its highly anisotropic feature even at a temperature where the
molecules have enough kinetic energy to become ‘mobile. It is then
possible that the translational topological defects, dislocations, have
lower excitation energy than the orientational topological defects, and
dislocations develop along with the orientational order. In this
high-pressure region, the dislocation-mediated melting theory for

anisotropic layers proposed by Ostlund and Halperin(4-28)

will probably
be applicable. The system may melt to an anisotropic liquid through a
continuous transition as predicted by this theory, then at a higher
temperature the system undergoes another transition to become an
isotropic liquid. Even if the orientational topological defects play an
important role in above transitions, due to the higher vortex
core-energy caused by smaller separation between molecules, the

transitions are possibly still continuous. If it is the true situation,

then we have a model liquid crystal system with the well-known

111

Lennard-Jones potential and an external parameter, namely pressure which
can be easily adjusted in computer simulation to obtain this phase.
Such a model should be ideal in the theoretical study of the properties
of 2D liquid crystal phase. From the view point of phase transition,
another interesting problem is then what happens between these zero
pressure and high pressure extreme cases? Multicritical points are
possibly present where phase boundary associated with the first order

transition meets critical line.

112

8.3.00 2,“ a"... 1‘.“ up.» new 15.00 13.50 29.00 g.“

 

 

 

 

 

V

3.

0" -6

Reduced Taper-am Is 0.?

8 a
3.1 Step: e-e From am To 1m .2
:1 hating From T - 0.68
% a
E“ 2:
z
5?.

'53. 8.
a“ 'u
2'
:3
a 3
g
o

8

o" "o.

3 8

fimo 2'.“ 5'.“ f.“ ban h. . . :

1 R 1 so than 1‘1.“ 25 u a u 6%
Fig. 4.12 The radial distribution function of the mass

centers of molecules at T=0.70, averaged over 3000 time
steps (60 ps). It shows liquid-like feature.

113

 

 

Fig. 4.13 A snapshot of the configuration of molecules at
T=0.70. The time step is 12000. There are many
dislocations.

(1.1)

(1.2)

(1.3)

(1.4)

(1.5)

(1.6)

(1.7)

(1.8)

(1.9)

(1.10)

(1.11)

(1.12)

(1.13)

114

REFERENCES FOR CHAPTER 1

One may gain an appreciation for the multitude of novel 2D
phenomena and the excitements in this field by referring to
Procedings of the International Conference on Ordering in Two
Dimensions, Lake Geneva, may 1980, Edited by S. R. Sinha
(North-Holland, Amsterdam, 1980)

J. Needham, Science and Civilization in China, Vol. 3, (Cambridge
University Press, 1959).

Phase Transitions in Surface Film, edited by J. G. Dash and
J. P. Ruvalds, (Plenum Press, NY 1979).

Recent experiments found that N monolayer was possiblly
incommensurate at high coverage. ee 5. C. Fain, Jr., and
H. You, Bull.Am.Phys.Soc. 19, 333 (1985).

For methane, R. Marx, and E. F. Wassermann, Surf.Sci. 111, 267
(1982); for ethane, P. Y. Hansen, R. Wang, H. Taub, H. Schechter,
D. G. Reichel, H. R. Danner, and G. P. Alldredge, Phys.Rev.Lett.
11, S72 (1984).

B. H. Grier, L. Passell, J. Eckert, H. Patterson, D. Richter, and
R. J. Rollefson, Phys.Rev.Lett. §§, 814 (1984); S. Nose, and
M. L. Klein, Phys.Rev.Lett. §§, 818 (1984).

J. P. MacTague, and M. Nielson, Phys.Rev.Lett. 11, 596, (1976);
M. Nielson, and J. P.MacTague, Phys.Rev. 212, 3096, (1979).

See, for example, R. Marx, and B. Christoffer, Phys.Rev.Lett.
11L790 (1983).

P. W. Stephens, P. A. Heiney, R. I. Birgeneau, P. M. Horn,
J. Stoltenbe and O. E. Vilches, Phys.Rev.Lett. £§11959 (1980).

R. D. Etters, P. Pan, and V. Chandrasekharan, Phys.Rev.Lett.
111645 (1980).

M. F. Toney, R. D. Diehl, and S. C. Fain, Phys.Rev. 827,6413
(1983).

J. M. Kosterlitz, and D. J. Thouless, J.Phys. 9;, L124, (1972);
J.Phys. §§, 1181, (1973).

B. I. Halperin, and D. R. Nelson, Phys.Rev.Lett. 11, 121 (1978);
D. R. Nelson, and B. I. Halperin, Phys.Rev. 112, 2457 (1979).

(1.14)

(1.15)

(1.16)
(1.17)

(1.18)

(1.19)

(1.20)

115

A. P. Young, Phys.Rev. 112, 1855 (1979).

Abraham showed that his computer simulation results suggested
first-order transition, 5. W. Koch, and F. F. Abraham, Phys.Rev.
9.2.7. 2964 (1983). A experimental example is, p. A. Heiney,
R. I. Birgeneau, C. .G.8rown, P. M. Horn, D. E. Moncton, and
P. W. Stephens, Phys.Rev.Lett. I1§L104 (1982).

S. Ostlund, and B. I. Halperin, Phys.Rev. 111, 335 (1981).
Y. Saito, Phys.Rev. 219, 6239 (1982).

H. Margenau, and N. R. Kestner, Theory of Intermolecular Forces,
(Pergamon Press, NY 1971).

A. I. Kitaigorodsky, Molecular Crystal and molecules, (Academic
Press, New YorK, 1973).

C. A. English, J. A. Venables, and D. R. Salahub,
Proc.Roy.Soc.London, A340, 81 (1974).

(2.1)

(2.2)

(2.3)

(2.4)

(2.5)

(2.6)

(2.7)

(2.8)

(2.9)

(2.10)

116

REFERENCES FOR CHAPTER 2

D. A. Goodings and M. Henkelman, Can.J.Phys. 12, 2898 (1971).

M. P. Toney, R. D. Diehl, and S. C. Pain, Phys.Rev. 211, 6413
(1983); P. A. Heiney, P. W. Stephens, S. G. J. Mochrie,
J. Akimatsu, and R. J. Birgeneau, Surf.Sci. 122, 539 (1983); and
the chapter 4 of this thesis.

R. K. Kalia, P. Vashishta, and S. D. Mahanti, Phys.Rev.Lett. 12,
676 (1982).

For parameters of interest in this chapter, we find that the
next-nearest neighbors contribute less than 10% of the
nearest-neighbor contribution.

J. M. Luttinger and L. Tisza, Phys.Rev. 12, 954 (1946);
J. M. Luttinger, ibid. 21, 1015 (1951).

The solution to the weak constraint problem satisfies the strong
constraint, thereby the exact ground state of the Hamiltonian.

See H. E. Stanley, Phase Transitions and Critical Phenomena,
(Oxford University Press, London, 1971).

P. Wegner, z.Phys. 222, 465 (1967).

For a system with Duo and K<O, J<O, and K/Jxo, we find that
n~exp[-const(T/J)|ln(K/J)|].

The actual Hamiltonian, when terms of order a6 are kept, has a
form much more complicated than Eq. [7]. However, the additional
terms drop out in H beyond those given in Eq. [7] for ordered
configurations shoag in Fig. 2.2. Therefore dropping these
additional terms in H and studying changes in parameters D, J
and x is meaningful Bfily when one wants to test thg stability of
the ordered ground-state configurations when O(a ) terms are
included.

(3.1)

(3.2)

(3.3)

(3.4)

(3.5)

(3.6)

(3.7)

(3.8)

(3.9).

(3.10)

(3.11)

(3.12)

(3.13)

117

REFERENCES FOR CHAPTER 3

Procedings of the International Conference on Ordering in Two
Dimensions, Lake Geneva, may 1980, Edited by S. K. Sinha
(North-Holland, Amsterdam, 1980)

When the substrate potential is strong and the density of
adsorbed molecules is not too large, the principle axes of the
molecules orient parallel to the surface; 91 is the orientation
of the i-th molecule.

The previous chapter of this thesis.

0. G. Mouritsen, and A. J. Berlinsky, Phys.Rev.Lett. 12, 181
(1982). K. Kaski, S. Knmar, J. D. Gunton, and P. A. Rikvold,
Phys.Rev. 222, 4420 (1984).

N. D. Mermin, and H. Wagner, Phys.Rev.Lett. 11, 1133 (1966).

H. E. Stanley, and T. Kaplan, Phys.Rev.Lett. 11, 913 (1966);
H. E. Stanley, Phys.Rev.Lett. 22, 150 (1968).

J. M. Kosterlitz, and D. J. Thouless, J.Phys. 22, 1181 (1973);
J. M. Kosterlitz, J.Phys. 21, 1046 (1974).

P. Wagner, z.Phys. 222, 465 (1967).

V. Jose, L. P. Radanoff, S. Kirkpatrick, and D. R. Nelson,
Phys.Rev. 212, 1217 (1977).

M. B. Einhorn, R. Savit, and E. Ravinovici, Nucl.Phys. 112, 16
(1980); R. Savit, Rev.Mod.Phys. 22, 453 (1980).

J.Tobochnik, and G. V. Chester, Phys.Rev. 222, 3761 (1979).

G. C. Grest, and D. J. Srolovitz, Phys.Rev. 222, 2693 (1983); In
addition T. Schneider and E. Stoll have found that for 2D HY
model with quartic anisotropy vortices were adsent near Tc' see
Solitons and Condensed Matter Physics, (Springer-Verlag, 1978,
p. 154).

A. A. Migdal, Sov.Phys. JETP 12, 743 (1976); L. P. Kadanoff,
Ann.Pyhs. (N.Y.) 122, 359 (1976).

118

(3.14) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Zelle,
and E. Zeller, J.Chem.Phys. 21, 1089 (1953).

(3.15) S. H. Schenker, and J. Tobochnik, Phys.Rev. 222, 4462 (1980).
(3.16) N. D. Mermin, Rev.Mod.Phys. 21, 591 (1971).

(3.17) B. M. McCoy, and T. T. Wu, The Two-Dimensional Ising Model,
(Harvard University Press, Cambridge, Massachusetts 1973).

(3.18) C. Domb, in Phase Transitions and Critical Phenomena,
Vol. 3,(Academic Press, 1975).

(4.1)

(4.2)

(4.3)

(4.4)

(4.5)

(4.6)

(4.7)

(4.8)

(4.9)

(4.10)

(4.11)

(4.12)

(4.13)

(4.14)

119

REFERENCES FOR CHAPTER 4

F. Schwabl, Perroelectrics, 21, 171 (1980); S. D.Mahanti,
D. Sahu, Phys.Rev.Lett. 12, 936 (1982).

M. P.Toney, R. D. Diehl, and S. C. Pain, Phys.Rev. 221, 6413
(1983); P. A. Heiney, P. W. Stephens, S. G. J. Mochrie,
J. Akikatsu, and R. J. Birgeneau, Surf.Sci. 122, 539 (1983).

R. D. Diehl, and S. C. Fain, Surf.Sci. 122, 116 (1983).

Ordering in Strongly Fluctuating Condensed Matter Systems, edited
by T. Riste, (Plenum Press, NY 1979).

For example, R. Max, and B. Christoffer, Phys.Rev.Lett. 2_
(1983).

P
Q

790
O. G. Mouritsen, and A. J. Berlinsky, Phys.Rev.Lett. 12, 181
(1982).

R. D. Etters, R. Pan, and V. Chandrasekharan, Phys.Rev.Lett. 12,
645 (1980); R.Pan, R. D. Etters, K. Nobashi, and
V. Chandrasekharan, J.Chem.Phys. 11, 1035 (1982).

See Chapter 2.

E. Domany, and E. K. Riedel, Phys.Rev.Lett. I 12, 561 (1978);
Phys.Rev. 212, 5817 (1979).

R. Savit, Rev.Mod.Phys. 22, 453 (1980).

J. M. Kosterlitz, and D. J. Thouless, J.Phys. 22, 1181 (1973);

'J. M. Kosterlitz, J.Phys. 91, 1046 (1974).

Monte Carlo Methods in Statistical Physics , Edited by K. Binder
(Springer-Verlag, Berlin, 1978); Applications of the Monte Carlo
Method in Statistical Physics , Edited by R. Binder
(Springer-Verlag, Berlin, 1980).

H. C. Anderson, J.Chem.Phys. 12, 2384 (1980).

M. Perrinello, and A. Rahman, Phys.Rev.Lett. 12, 1196 (1980);
J.App1.Phys. 21, 7182, (1980); J.Chem.Phys. lg, 2662 (1980).

(4.15)

(4.16)

(4.17)

(4.18)

(4.19)

(4.20)

(4.21)

(4.22)

(4.23)
(4.24)

(4.25)

(4.26)

(4.27)

(4.28)

(4.29)

120

S. Nose, and M. L. Klein, J.Chem.Phys. 12, 6928 (1983);
Mo1.Phys. 22, 1055 (1983).

C. W. Gear, Numerical Initial Value Problems in Ordinary
Differential Equations, (Pretice Hall, Englewood Cliffs, N.J.
1971).

F. James, Rep.Prog.Phys. 12, 73 (1980).

See, for example, L. Verlet, Phys.Rev. 122, 98 (1967).

M. Sprik, R. W. Impey, and M. L. Klein, Phys.Rev. 222, 4368

(1984);

M. Born, and K. Huang, Dynamical Theory of Crystal Lattice,
(Oxford University Press, London, 1968)
G. C. DeFotic, Phys.Rev. B23, 4714 (1981).

M. B. Einhorn,
(1980).

R. Savit, and E. Ravinovici, Nuc1.Phys. 112, 16

Y. Saito, Phys.Rev. 22, 6239 (1982).

S. W. Koch, and P. P. Abraham, Phys.Rev. 221, 2964 (1983).
R. K. Kalia, P. Vashishta, and S. D. Mahanti, Phys.Rev.Lett. 12,
676 (1982).

J. E. Van Himbergen, Phys.Rev.Lett. 22, 5 (1984). However,
there is other opinion about the nature of the transition in this
model, see H. J. F. Knops, Phys.Rev. 222, 470 (1984).

v. Jose, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson,
Phys.Rev. 212, 1217 (1977).

See, for example, R. Swendsen, Phys.Rev.Lett. 12, 1302 (1982);
E. Domany, M. Schick, and R. H. Swenden, Phys.Rev.Lett. 22, 1535

(1984).

S. Ostlund, and B. I. Halperin, Phys.Rev. 222, 335 (1981).

"mmmm‘s