Spin dynamics simulations of collective excitations in magnetic liquids
A novel approach is developed for computer simulation studies of dynamical properties of spin liquids. It is based on the Liouville operator formalism of Hamiltonian dynamics in conjunction with Suzuki-Trotter-like decompositions of exponential propagators. As a result, a whole set of symplectic t...
Збережено в:
| Дата: | 2000 |
|---|---|
| Автори: | , , |
| Формат: | Стаття |
| Мова: | English |
| Опубліковано: |
Інститут фізики конденсованих систем НАН України
2000
|
| Назва видання: | Condensed Matter Physics |
| Онлайн доступ: | https://nasplib.isofts.kiev.ua/handle/123456789/121006 |
| Теги: |
Додати тег
Немає тегів, Будьте першим, хто поставить тег для цього запису!
|
| Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Цитувати: | Spin dynamics simulations of collective excitations in magnetic liquids / I.P. Omelyan, I.M. Mryglod, R. Folk // Condensed Matter Physics. — 2000. — Т. 3, № 3(23). — С. 497-514. — Бібліогр.: 29 назв. — англ. |
Репозитарії
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
nasplib_isofts_kiev_ua-123456789-121006 |
|---|---|
| record_format |
dspace |
| spelling |
nasplib_isofts_kiev_ua-123456789-1210062025-02-09T20:49:11Z Spin dynamics simulations of collective excitations in magnetic liquids Спiн-динамiчнi розрахунки колективних збуджень у магнiтних piдинах Omelyan, I.P. Mryglod, I.M. Folk, R. A novel approach is developed for computer simulation studies of dynamical properties of spin liquids. It is based on the Liouville operator formalism of Hamiltonian dynamics in conjunction with Suzuki-Trotter-like decompositions of exponential propagators. As a result, a whole set of symplectic time-reversible algorithms has been introduced for numerical integration of the equations of motion at the presence of both translational and spin degrees of freedom. It is shown that these algorithms can be used in actual simulations with much larger time steps than those inherent in standard predictor-corrector schemes. This has allowed one to perform direct quantitative measurements for spin-spin, spin-density and density-density dynamical structure factors of a Heisenberg ferrofluid model for the first time. It was established that like pure liquids the density spectrum can be expressed in terms of heat and sound modes, whereas like spin lattices in the ferromagnetic phase there exists one primary spin in the shape of spin- spin dynamic structure factors describing the longitudinal and transverse spin fluctuations. As it was predicted in our previous paper [Mryglod I., Folk R. et al., Physica A277 (2000) 389] we found also that a secondary wave peak appears additionally in the longitudinal spin-spin dynamic structure factor. The frequency position of this peak coincides entirely with that for a sound mode reflecting the effect of the liquid subsystem on spin dynamics. The possibility of longitudinal spin wave propagation in magnetic liquids at sound frequency can be considered as a new effect which has yet to be tested experimentally. Пpопонується новий пiдхiд до дослiдження динамiчних властивостей спiнових piдин у комп’ютеpному експеpиментi. Пiдхiд базується на фоpмалiзмi опеpатоpiв Лiувiлля для опису Гамiльтонової динамiки у поєднаннi з методом фактоpизацiї експонентних пpопагатоpiв типу Сузукi-Тpоттеpа. У pезультатi розвинуто сукупність нових алгоритмів для числового iнтегpування piвнянь pуху за наявностi як тpасляцiйних, так i спiнових ступеней вiльностi, які зберігають фазовий об’єм і є часо-зворотніми. Показано, що цi алгоpитми можуть викоpистовуватися у комп’ютеpних розрахунках із набагато бiльшими часовими кpоками порівняно з тими, якi пpитаманнi стандаpтним пpогнозокоpектованим схемам. Це нам дало змогу впеpше виконати прямі розрахунки динамiчних стpуктуpних фактоpiв типу спiн-спiн, спiн- густина i густина-густина для Гайзенбеpгiвської моделi феpофлюїду. Встановлено, що подiбно до чистих piдин густинний спектp може бути описаний в термінах теплової та звукових мод. Водночас, подiбно до спінових гpаткових систем,у феромагнітній фазі у поздовжньому та попеpечному динамічних структурних факторах спін-спін як основний спостерігається спін-хвильовий пік. Окрім того, у поздовжньому динамічному структурному факторі спін-спін виявлено додатковий пік, поява якого пеpедбачалася нами раніше [Mryglod I., Folk R. et al., Physica A277 (2000) 389]. Частотна позиція цього піку повністю збігається з частотою звукових збуджень, що відображає вплив рідинної підсистеми на спінову динаміку. Можливiсть пошиpення поздовжнiх спiнових хвиль на звуковій частоті слід розглядати як новий ефект, що потребує експериментальної перевірки. Part of this work was supported by the Fonds zur Förderung der wissenschaftlichen Forschung under Project No. P12422–TPH. 2000 Article Spin dynamics simulations of collective excitations in magnetic liquids / I.P. Omelyan, I.M. Mryglod, R. Folk // Condensed Matter Physics. — 2000. — Т. 3, № 3(23). — С. 497-514. — Бібліогр.: 29 назв. — англ. 1607-324X DOI:10.5488/CMP.3.3.497 PACS: 02.60.Cb, 75.40.Gb, 75.50.Mm, 76.50.+g, 95.75.Pq https://nasplib.isofts.kiev.ua/handle/123456789/121006 en Condensed Matter Physics application/pdf Інститут фізики конденсованих систем НАН України |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| language |
English |
| description |
A novel approach is developed for computer simulation studies of dynamical properties of spin liquids. It is based on the Liouville operator formalism
of Hamiltonian dynamics in conjunction with Suzuki-Trotter-like decompositions of exponential propagators. As a result, a whole set of symplectic
time-reversible algorithms has been introduced for numerical integration of
the equations of motion at the presence of both translational and spin degrees of freedom. It is shown that these algorithms can be used in actual
simulations with much larger time steps than those inherent in standard
predictor-corrector schemes. This has allowed one to perform direct quantitative measurements for spin-spin, spin-density and density-density dynamical structure factors of a Heisenberg ferrofluid model for the first time.
It was established that like pure liquids the density spectrum can be expressed in terms of heat and sound modes, whereas like spin lattices in
the ferromagnetic phase there exists one primary spin in the shape of spin-
spin dynamic structure factors describing the longitudinal and transverse
spin fluctuations. As it was predicted in our previous paper [Mryglod I.,
Folk R. et al., Physica A277 (2000) 389] we found also that a secondary
wave peak appears additionally in the longitudinal spin-spin dynamic structure factor. The frequency position of this peak coincides entirely with that
for a sound mode reflecting the effect of the liquid subsystem on spin dynamics. The possibility of longitudinal spin wave propagation in magnetic
liquids at sound frequency can be considered as a new effect which has
yet to be tested experimentally. |
| format |
Article |
| author |
Omelyan, I.P. Mryglod, I.M. Folk, R. |
| spellingShingle |
Omelyan, I.P. Mryglod, I.M. Folk, R. Spin dynamics simulations of collective excitations in magnetic liquids Condensed Matter Physics |
| author_facet |
Omelyan, I.P. Mryglod, I.M. Folk, R. |
| author_sort |
Omelyan, I.P. |
| title |
Spin dynamics simulations of collective excitations in magnetic liquids |
| title_short |
Spin dynamics simulations of collective excitations in magnetic liquids |
| title_full |
Spin dynamics simulations of collective excitations in magnetic liquids |
| title_fullStr |
Spin dynamics simulations of collective excitations in magnetic liquids |
| title_full_unstemmed |
Spin dynamics simulations of collective excitations in magnetic liquids |
| title_sort |
spin dynamics simulations of collective excitations in magnetic liquids |
| publisher |
Інститут фізики конденсованих систем НАН України |
| publishDate |
2000 |
| url |
https://nasplib.isofts.kiev.ua/handle/123456789/121006 |
| citation_txt |
Spin dynamics simulations of collective excitations in magnetic liquids / I.P. Omelyan, I.M. Mryglod, R. Folk // Condensed Matter Physics. — 2000. — Т. 3, № 3(23). — С. 497-514. — Бібліогр.: 29 назв. — англ. |
| series |
Condensed Matter Physics |
| work_keys_str_mv |
AT omelyanip spindynamicssimulationsofcollectiveexcitationsinmagneticliquids AT mryglodim spindynamicssimulationsofcollectiveexcitationsinmagneticliquids AT folkr spindynamicssimulationsofcollectiveexcitationsinmagneticliquids AT omelyanip spindinamičnirozrahunkikolektivnihzbudženʹumagnitnihpidinah AT mryglodim spindinamičnirozrahunkikolektivnihzbudženʹumagnitnihpidinah AT folkr spindinamičnirozrahunkikolektivnihzbudženʹumagnitnihpidinah |
| first_indexed |
2025-11-30T16:04:10Z |
| last_indexed |
2025-11-30T16:04:10Z |
| _version_ |
1850231911390117888 |
| fulltext |
Condensed Matter Physics, 2000, Vol. 3, No. 3(23), pp. 497–514
Spin dynamics simulations of collective
excitations in magnetic liquids
I.P.Omelyan 1 , I.M.Mryglod 1,2 , R.Folk 2
1 Institute for Condensed Matter Physics
of the National Academy of Sciences of Ukraine,
1 Svientsitskii Str., 79011 Lviv, Ukraine
2 Institute for Theoretical Physics, University of Linz, A-4040 Linz, Austria
Received April 29, 2000
A novel approach is developed for computer simulation studies of dynami-
cal properties of spin liquids. It is based on the Liouville operator formalism
of Hamiltonian dynamics in conjunction with Suzuki-Trotter-like decompo-
sitions of exponential propagators. As a result, a whole set of symplectic
time-reversible algorithms has been introduced for numerical integration of
the equations of motion at the presence of both translational and spin de-
grees of freedom. It is shown that these algorithms can be used in actual
simulations with much larger time steps than those inherent in standard
predictor-corrector schemes. This has allowed one to perform direct quan-
titative measurements for spin-spin, spin-density and density-density dy-
namical structure factors of a Heisenberg ferrofluid model for the first time.
It was established that like pure liquids the density spectrum can be ex-
pressed in terms of heat and sound modes, whereas like spin lattices in
the ferromagnetic phase there exists one primary spin in the shape of spin-
spin dynamic structure factors describing the longitudinal and transverse
spin fluctuations. As it was predicted in our previous paper [Mryglod I.,
Folk R. et al., Physica A277 (2000) 389] we found also that a secondary
wave peak appears additionally in the longitudinal spin-spin dynamic struc-
ture factor. The frequency position of this peak coincides entirely with that
for a sound mode reflecting the effect of the liquid subsystem on spin dy-
namics. The possibility of longitudinal spin wave propagation in magnetic
liquids at sound frequency can be considered as a new effect which has
yet to be tested experimentally.
Key words: Hejsenberg fluids, spin relaxation, collective modes,
symplectic integration, Suzuki-Trotter decompositions
PACS: 02.60.Cb, 75.40.Gb, 75.50.Mm, 76.50.+g, 95.75.Pq
1. Introduction
Magnetic systems have been the subject of numerous theoretical and computer
experimental studies [1–6]. These studies were devoted mainly to describe static
properties, such as phase transitions, scaling, and critical phenomena. The theoret-
ical description of dynamic properties presents a more difficult problem which until
c© I.P.Omelyan, I.M.Mryglod, R.Folk 497
I.P.Omelyan, I.M.Mryglod, R.Folk
now cannot be solved quantitatively even for the simplest magnetic models. For
this reason, the method of molecular dynamics (MD) was used over recent years to
investigate dynamic critical behaviour, dynamic scaling, and processes of spin relax-
ation. However, the MD simulations were restricted to lattice systems exclusively,
in particular, to classical XY and Heisenberg models in d = 2 as well as Heisenberg
ferro- and antiferromagnets in d = 3 dimensions [7–11].
The need in extensions of MD simulations to continuum disordered systems is
motivated by a great variety of additional physical phenomena arising when both
spin (orientational) and liquid (translational) degrees of freedom are taken into
account. The existing theoretical investigations of magnetic liquid dynamics dealt to
a great extent with phenomenological approaches [12–14]. Recently, the equations
for time correlation functions and collective mode spectra of a Heisenberg ferrofluid
model have been consequently derived using the method of nonequilibrium statistical
operator [15,16]. As a result, a simple model of spin relaxation was considered, and
the frequency matrix as well as the matrix of memory functions were calculated.
However, this approach itself requires some additional information on the dynamic
behaviour expressed in terms of correlation times. It is also worth mentioning that
the main part of the theoretical results was obtained for a spin subsystem only and
the mutual effect of the liquid subsystem on the spin dynamics has not been studied
in detail.
In view of the afore said, computer experiments should be considered as the chief
tool for quantitative studies of the dynamic properties. Nevertheless, there were no
attempts known to simulate magnetic liquids within the MD approach. This can
be explained by the absence of a suitable algorithm for solving the corresponding
equations of motion. The traditional numerical methods [17] cannot be used for
observations over a many-body system because of their high instability on MD scales
of time. It was established for lattice systems that even standard predictor-corrector
schemes appear to be inefficient because of very poor conservations of energy and
spin lengths [18].
Obviously, an efficient algorithm for long-duration integrations should be stable,
accurate, fast and easy to implement. There is only a limited group of integrators
satisfying these requirements. Among them one can distinguish the velocity Verlet
algorithm [19,20] which allows us to obtain a high level of accuracy with mini-
mal computational costs measured in terms of time-consuming force evaluations.
However, the velocity Verlet and other similar schemes [20–22] were constructed to
integrate translational motion without the presence of any orientational degrees of
freedom. In our case the pattern is still more complicated since the translational
positions and momenta are coupled with spin orientations in a characteristic way
and, therefore, all these dynamical variables should be integrated simultaneously.
This requires substantial revision of the liquid dynamic algorithms.
Quite recently, new spin dynamics algorithms have been proposed to simulate
lattice systems [18]. They were constructed within the Suzuki-Trotter technique [23]
using sublattice decompositions of spin dynamics. As was demonstrated, the sublat-
tice decomposition algorithms unlike usual schemes exactly conserve individual spin
498
Spin dynamics simulations. . .
lengths, exhibit good stability and energy conserving properties, and allow much
larger time steps than predictor-corrector methods. The main problem with these
algorithms lies in the fact that they are applicable only to spin systems when the
decomposition on two or several noninteracting sublattices is possible. Such algo-
rithms cannot be used for models with arbitrary spatial spin distributions and, thus,
not for spin liquids.
In the current paper we develop the Suzuki-Trotter decomposition method and
derive the desired spin liquid algorithms of different orders in the time step. The
algorithms are tested in actual MD simulations on a Heisenberg ferrofluid model and
compared with predictor-corrector schemes. The computations of spin-spin, spin-
density and density-density time correlation functions are also carried out and the
discussion on the surprising results obtained is added.
2. Basic equations of motion
Let us consider a classical fluid composed of N magnetic particles of massm each
containing continuous three-component spin si = (sxi , s
y
i , s
z
i ) with the fixed length
|si| = 1 for each site i = 1, 2, . . . , N . A typical model Hamiltonian for such a system
can be presented as
H =
N
∑
i=1
mvi
2
2
+
N
∑
i<j
(
ϕij − Jij(s
x
i s
x
j + syi s
y
j + λszi s
z
j )
)
−G
N
∑
i=1
szi . (1)
Here vi ≡ dri/dt and ri are the translational velocity and position of particle i, re-
spectively, ϕij ≡ ϕ(rij) is the liquid potential, Jij ≡ J(rij) denotes the exchange
integral corresponding to a pair (i, j) of spins with the interparticle separation
rij = |ri − rj|, the exchange anisotropy parameter is denoted by λ, and G is the in-
tensity of an external spatially and timely homogeneous magnetic fieldG ≡ (0, 0, G)
directed along the z-th axis. Equation (1) represents, in fact, the isotropic (λ = 1)
or anisotropic (λ 6= 1) Heisenberg ferro- or the corresponding antiferro-fluid for
Jij > 0 and Jij < 0, respectively. Although the main results which will be obtained
in this paper are applicable for a larger class of Hamiltonians, we restrict ourselves
to simplify notations to the basic model (1).
In the Liouville formulation of Hamiltonian dynamics, the equations of motion
for the set ρ ≡ {ri,vi, si} of microscopic variables are of the form
dρ
dt
= [ρ, H ] ≡ Lρ(t) , (2)
where
L =
N
∑
i=1
(
vi·
∂
∂ri
+
fi
m
·
∂
∂vi
− 1
~
[
(gi +G)×si
]
·
∂
∂si
)
≡ Lr + Lv + Ls (3)
is the Liouville operator [15] of the system with
fi = −
N
∑
j(j 6=i)
[dϕij
drij
− dJij
drij
(
sxi s
x
j + syi s
y
j + λszi s
z
j
)]rij
rij
(4)
499
I.P.Omelyan, I.M.Mryglod, R.Folk
being the translational force, and
gi ≡ (gxi , g
y
i , g
z
i ) =
N
∑
j(j 6=i)
Jij(s
x
j , s
y
j , λs
z
j) (5)
the internal magnetic field. Note that the components Lr, Lv and Ls act only on
position, velocity and spin, respectively, and the quantum Poisson bracket [ , ] was
used to derive the expression for Ls. If an initial configuration ρ(0) is determined,
the time evolution of ρ(t) can be obtained numerically by integrating equation (2).
Taking into account the symmetries ϕij = ϕji and Jij = Jji, it follows from
equations (1) and (2) that the total energy E ≡ H and the total momentum P =
m
∑
i vi of the system are integrals of motion, i.e. dE/dt = 0 and dP/dt = 0. These
symmetries also impose a conservation law for the magnetization M =
∑
i si of the
isotropic Heisenberg model at the absence of an external magnetic field. For the
anisotropic case (when λ 6= 1 and/or G 6= 0) only the z-th component Mz of M
will not change in time. The structure of equation (2) requires the conservation of
individual spin lengths as well. Moreover, it can be shown readily that our equations
of motion are time reversible, i.e., they are invariant with respect to the set t →
−t, {vi, si} → {−vi,−si} of transformations. Besides, the exact solutions behave
symplectically, i.e., the phase volume of {r i,vi} will remain constant. Obviously, all
the above properties cannot be perfectly reproduced during the integration within
any numerical scheme. This is a typical situation in MD simulations of liquids.
However, as we shall show in the next section, some of the integrals of motion
as well as the time-reversibility and symplecticity can be maintained although the
numerical trajectories are generated with a limited accuracy.
3. Algorithms for numerical integration
Within the Liouville exponential operator formalism, the solutions to the equa-
tions of motion (2) are cast in the form
ρ(t+ h) = eLhρ(t) = e(Lr+Lv+Ls)hρ(t) , (6)
where h denotes the time step. Since the exponential propagator eLh cannot be
evaluated exactly, it is necessary to introduce some approximations. These approx-
imations should take into account the relatively small time step h actually used in
MD simulations.
3.1. Second-order version
Assuming for the moment that spin variables are frozen, i.e. letting L s → 0,
one comes to the usual (liquid-like) equations of motion. They can be solved in a
quite efficient way applying the second-order VV integrator [19,20] which is based on
the Suzuki-Trotter (ST) formula e(Lr+Lv)h = eLvh/2eLrheLvh/2 + O(h3). Unfreezing
now the spin variables and treating the sum Lr + Ls as one operator we obtain
500
Spin dynamics simulations. . .
immediately: e(Lr+Lv+Ls)h = eLvh/2e(Lr+Ls)heLvh/2+O(h3) (note that the ST formula
is valid for arbitrary two operators). The spin-position subpropagator can further
be decomposed in a similar way, e(Lr+Ls)h = eLrh/2eLsheLrh/2 +O(h3), resulting into
the full propagation
ρ(t+ h) = eLv
h
2 eLr
h
2 eLsheLr
h
2 eLv
h
2ρ(t) +O(h3) ≡ D(t, h)ρ(t) +O(h3) . (7)
Of course, other decompositions are also possible, but then the local fields g i or/and
forces fi will need to be updated (the most time-consuming operations) more fre-
quently, thus reducing the efficiency of the computations.
The main idea of the decompositions consists in obtaining subpropagators which
can be evaluated analytically. As can be shown [20], the position eLrτ =
∏
i e
Lri
τ
and velocity eLvτ =
∏
i e
Lvi
τ propagation present shift operators, namely,
eLri
τri = ri + viτ , eLvi
τvi = vi +
fi
m
τ, (8)
where Lri
= vi·∂/∂ri, Lvi
= fi/m·∂/∂vi, and the equalities Lri
Lrj
= Lrj
Lri
,
Lvi
Lvj
= Lvj
Lvi
have been used. It is worth emphasizing that the components
Lr and Lv (as well as Ls) do not commute between themselves. For this reason, the
position and velocity shifts (8) must be performed in the rigorous order as specified
by equation (7) and applied to actual current values of ri and vi within the time
step.
The spin subdynamics is described in equation (7) by the exponential operator
eLsh, where Ls =
∑
i Lsi
and Lsi
= [ωi×si]·∂/∂si with ωi = −(gi +G)/~ being the
local Larmor frequency. This operator has no simple explicit form, given that the
internal magnetic field gi, acting on each particle, depends on the orientations of
generally speaking (see equation (5)) all other spins of the system. In other words,
it is itself time dependent and is, therefore, not known a priori. A possibility to
obtain the explicit solutions lies in the following. Since all the partial components
Lsi
do not commute each other, it is quite natural to find an ST-like decomposition
for the whole set (i = 1, . . . , N) of these noncommutative operators. Consecutively
(k = 1, . . . , N −1) dividing the subsets (i = k, . . . , N) into groups composed of only
two operators (Lsk
and
∑N
j=k+1Lsj
) and using the usual ST formula for them, one
obtains
eLsh = eLs1
h
2 . . . eLsN−1
h
2 eLsN
heLsN−1
h
2 . . . eLs1
h
2 +O(h3) . (9)
Equation (9) constitutes an ST analogue for arbitrary number of operators and
is accurate to the same order as that of the initially truncated (see equation (7))
terms O(h3). Again, other O(h3) decompositions can be introduced. However, only
equation (9) will lead to a scheme with minimal number of local field recalculations.
The problem is now much simplified because, according to equation (9), each
current value of si is updated spin by spin at a fixed instantaneous Larmor frequency
ωi, and this case allows analytical solutions. The result is
eLsi
τsi(t) =
(
I+Wi sin(ωiτ) +W2
i
(
1− cos(ωiτ)
))
si(t) ≡ Θi(t, τ)si(t) , (10)
501
I.P.Omelyan, I.M.Mryglod, R.Folk
where Θi(t, τ) denotes an orthonormal (ΘΘ+ = I) matrix which rotates around axis
ωi on angle ωiτ , and Wi = W(ωi) is a skewsymmetric matrix (Wαβ = −Wβα) with
the elements WXY = −ωZ/ω, WXZ = ωY /ω, WY Z = −ωX/ω. The trigonometric
functions can be avoided in view of the fact that the above decompositions are
only correct within O(h3) uncertainty. It is, therefore, sufficient to replace them by
rational counterparts, cos ξ = (1−ξ2/4)/(1+ξ2/4)+O(ξ3) and sin ξ = ξ/(1+ξ2/4)+
O(ξ3), which maintain the orthonormality of Θ. Then the spin rotation reduces to
eLsi
τsi(t) =
si(t) + [ωi×si(t)]τ + 1
2
(
ωi(ωi·si(t))− 1
2
(ωi·ωi)si(t)
)
τ 2
1 + 1
4
(ωiτ)
2 +O(τ 3) . (11)
This completes our spin liquid dynamics (SLD) algorithm of the second-order.
It can be shown readily that the algorithm derived generates solutions which are
time reversible and symplectic. Indeed, the initial propagator was decomposed (see
equations (6), (7) and (9)) into subparts symmetrically, and, as a consequence, the
final expressions for ri, vi and si will be invariant with respect to the transformation
h → −h. Furthermore, simple shifts (8) applied separately in position and velocity
space do not change the phase volume. These properties are very important for
our purpose because, as is now well established [20,21], the stability of an algorithm
normally follows from its time reversibility and symplecticity. Another nice property
of the SLD integrator is its exact conservation of individual spin lengths (rotations
given by equation (10) or (11) do not change the norm of vectors). This is crucial
as well, since the condition |si| = const represents a major part of the definition for
the model.
3.2. Higher-order versions
In many MD applications the accuracy achieved by second-order schemes is quite
sufficient to obtain reliable results. However, if very high precision is required, higher-
order algorithms may appear to be more efficient. The reason of this is that within
these algorithms the truncation errors decrease more rapidly with decreasing of the
step size. Applying very accurate MD integrators is especially important at phase
transitions, for long-time correlated properties or when investigating derivatives of
the thermodynamic functions.
A possibility to construct higher-order versions of the SLD algorithm lies in
employing a so-called multiple staging technique [23]. This technique was used earlier
for lattice spin dynamics [18] in the framework of usual ST decompositions. It was
realized by us that the multiple staging integration can be performed with slight
modifications within our approach as well. As a result, we have obtained the following
expression
ρ(t+ h) =
P
∏
p=1
D(t, ξph)ρ(t) +O(hq+1) (12)
for the q-th order propagation, where the coefficients ξp are chosen at a given number
P to provide the highest possible value of q. For example, the fourth-order (q = 4)
502
Spin dynamics simulations. . .
algorithm (SLD4) can be directly derived from equation (12) using P = 5 and
the coefficients ξ1 = ξ2 = ξ4 = ξ5 ≡ ξ = 1/(4 − 41/3), and ξ3 = 1 − 4ξ. That is
very interesting, these coefficients coincide with those found within the usual ST
decomposition of exponential operators [23]. At P = 1 and ξ1 = 1, equation (12)
reduces to the second order SLD algorithm.
It is worth remarking that within the higher-order versions (q > 2), the individual
spin propagation must be carried out using exact rotations (10) rather than their
second-order counterparts (11). Alternatively, the trigonometric functions can be
avoided in equation (10) again, but then more precise rational approximations should
be considered for them. In particular, the forth-order approximations are:
sin ξ =
ξ − 1
12
ξ3
1 + 1
12
ξ2 + 1
144
ξ4
+O(ξ5) , cos ξ =
1− 5
12
ξ2 + 1
144
ξ4
1 + 1
12
ξ2 + 1
144
ξ4
+O(ξ5) . (13)
Expressions (13) ensure the orthonormality of Θ and may be used within the SLD4
integration. Clearly, the fourth-order version SLD4 is time reversible since the mul-
tipliers ξp enter symmetrically into equation (12). Moreover, it conserves the phase
volume and individual spin lengths because the conservations are observed at each
stage p. The corresponding versions with q > 4 can also be explicitly introduced but
they appear to be too demanding for practical calculations.
3.3. Multiple time scale integration
The existence of fast and slow processes in the system allows one to improve
the efficiency of the decomposition integration using a multiple time scale (MTS)
approach. The main idea of this approach is based on using a short time step for the
rapidly varying strong forces, and a large time step for the smooth weak forces arising
from the interactions at long interparticle distances. Then the most time-consuming
slow forces will be recalculated less frequently saving significantly the computer
time. Note that the MTS technique was originally constructed [22] to integrate pure
liquid dynamics. Here we consider a question of how the MTS integration can be
adapted to our SLD/SLD4 algorithms.
A way to realize this consists in splitting the liquid potential ϕ ij into fast
(rij < R) and slow (rij > R) varying subparts which correspond to the repulsion
(dϕij/drij < 0) and attraction (dϕij/drij > 0) intervals of the interparticle distance
rij. Then the total translational force (4) can be presented as the sum fi = f
(f)
i + f
(s)
i
of two terms,
f
(f)
i = −
N
∑
j(j 6=i,rij<R)
dϕij
drij
rij
rij
,
f
(s)
i = −
N
∑
j(j 6=i,rij>R)
dϕij
drij
rij
rij
+
N
∑
j(j 6=i)
[dJij
drij
(
sxi s
x
j + syi s
y
j + λszi s
z
j
)]rij
rij
, (14)
associated with the fast f
(f)
i and slow f
(s)
i varying subforces, respectively. This results
into the separation Lvi
= L
(f)
vi + L
(s)
vi of the velocity part of the Liouville operator,
503
I.P.Omelyan, I.M.Mryglod, R.Folk
where L
(f)
vi = 1
m
f
(f)
i ·
∂
∂vi
and L
(s)
vi = 1
m
f
(s)
i ·
∂
∂vi
. Acting similarly as in the case of the
derivation of ST-like decompositions (7) and (9), taking into account the above sep-
aration of Lvi
, and collecting the fastest exponential operators inside the innermost
cycles we obtain
ρ(t+ h) = eL
(s)
v
h
2
[
eL
(f)
v
h
4k eLr
h
2k eL
(f)
v
h
4k
]k
eLsh
[
eL
(f)
v
h
4k eLr
h
2k eL
(f)
v
h
4k
]k
eL
(s)
v
h
2ρ(t) ≡ Dmts
ρ(t) ,
(15)
where eL
(s)
v τ =
∏
i e
L
(s)
vi
τ and eL
(f)
v τ =
∏
i e
L
(f)
vi
τ . During such a MTS integration the
high energy collisions between particles at rij 6 R are described more accurately by
decreasing the time step to h/k, where k > 1 is a parameter of the method.
Of course, the MTS integration requires additional computational efforts for
performing 2k recalculations of f
(f)
i . But they will be compensated completely by
using larger time steps h being within the same accuracy. On the other hand, using
the same time step for both the MTS (15) and usual (7) integrations, the precision of
the calculations can be improved considerably with very little computational costs.
This is so because the processor time needed to update f
(f)
i appears, as a rule, to
be much smaller (at appropriately chosen values for the multistep parameter k and
R) than that to evaluate long-ranged term f
(s)
i . Note that setting R → 0 (then
f
(f)
i → 0 and f
(s)
i ≡ fi), the MTS spin dynamics propagation (15) reduces to the
usual single time step SLD integration (7). It is also possible to implement the MTS
technique within the SLD4 integration. Therefore it necessary to formally replace in
equation (12) the propagator D(t, h) by their MTS analogue Dmts(t, h).
4. Applications to a Heisenberg ferrofluid model
Our MD studies were based on an isotropic (λ = 1) Heisenberg ferrofluid model
at the absence of an external magnetic field (G = 0). The Yukawa function [6] has
been used to describe the spin-spin interactions. This function was truncated at
Rc = 2.5σ and shifted to be zero at the truncation point to avoid force singularities,
J(r) =
ǫ
[
σ
r
exp
(
σ − r
σ
)
− σ
Rc
exp
(
σ − Rc
σ
)]
, r < Rc
0 , r > Rc
. (16)
The liquid subsystem was modelled by a soft-core potential [21],
ϕ(r) =
{
4ε[(σ/r)12 − (σ/r)6] + ε , r < 21/6σ ≡ R
0 , r > R
. (17)
Here σ denotes the diameter of particles, whereas ǫ and ε are the intensities of the
spin-spin and core-core interactions, respectively. The simulations were carried out
in the cubic box (employing periodic boundary conditions) with N = 1000 spins at
a reduced density n∗ = N
V
σ3 = 0.6, a reduced temperature T ∗ = kBT/ǫ = 1.5 <
T ∗
c (where T ∗
c ≈ 2.06 is the critical temperature of the system [24]), a reduced
504
Spin dynamics simulations. . .
core intensity ε/ǫ = 1, and a dynamical coupling parameter d = σ(mǫ)1/2/~ = 2.
This parameter presents, in fact, the ratio τtr/τs, where τtr = σ(m/ǫ)1/2 and τs =
~/ǫ are the characteristic time intervals of varying translational and spin variables,
respectively.
As far as we shall investigate a ferrophase and deal with a microcanonical
(NVES) ensemble, a non-zero magnetization of the system should be specified ad-
ditionally. This quantity was taken from our separate single MC simulation [24],
〈S〉0/N = 0.6536 ± 0.0001, where 〈 〉0 denotes the canonical averaging (the finite-
size effects were estimated to be less than 1%). The MC simulation was performed
using a biasing scheme [25] for sampling orientational degrees of freedom. Note
that such a magnetization could be determined in MD simulations too, but then
a canonical ensemble should be considered. This requires thermostatic forces and
virtual magnetic fields to be included in the equations of motion, but this is beyond
the scope of this paper.
All test runs were started from an identical well equilibrated configuration. The
recalculation of local magnetic fields during spin subdynamics propagation (9) took
approximately the same processor time as that of translational forces, spending in
total 0.5 sec per step on the Origin 2000 workstation. It is worth emphasizing that
contrary to pure spin dynamics [18] (when auxiliary MC cycles are involved to gen-
erate configurations as initial conditions for equations of motion), the equilibration
of our system can be performed within NVES MD simulations exclusively (where
the initial values for si should be generated at the very beginning of the equilibration
in order to fulfill the constraint S ≡ 〈S〉0). This is possible due to the presence of
energy exchange between the spin and liquid subsystems.
The MD results for the total energy E∗ = E/ǫ and total spin S as depending
on the length t of the simulations are presented in figure 1 (subsets (a)–(d) and
(e)–(h), respectively). Four time steps, namely, h∗ = h/τtr = 0.00125, 0.0025, 0.005
and 0.01, were used to integrate the equations of motion within our decomposition
SLD algorithm of the second order (equations (7)–(9) and (11)). For the purpose of
comparison, the results obtained by us using a well established Adams-Bashforth-
Moulton (ABM) predictor-corrector scheme [17] of the fourth-order are also included
in figure 1 (dashed curves in subsets (a) and (b)). We mention that in this scheme
the dynamical variables are first predicted, ρ(t + h) = ρ(t) + [55ρ̇(t) − 59ρ̇(t −
h) + 37ρ̇(t − 2h) − 9ρ̇(t − 3h)]h/24 + O(h5), and further iteratively corrected as
ρ(t+ h) = ρ(t) + [9ρ̇(t+ h) + 19ρ̇(t)− 5ρ̇(t− h) + ρ̇(t− 2h)]h/24 +O(h5), where
in our case ρ̇ = {vi, fi/m, [ωi×si]}. As can be seen from figure 1 (a), the ABM
integrator exhibits similar equivalence in energy conservation at the smallest time
step h∗ = 0.00125. However, this and other similar predictor-corrector schemes are
neither reversible nor symplectic, and there is no evidence that the extra order
obtained is relevant, since they have a very narrow region of stability and cannot be
used for larger step sizes. This is demonstrated in figure 1 (b) for the ABM integrator
(see dashed curve), where the total energy deviates systematically from the initial
level leading to the fluctuations at the end of the run which exceed the fluctuations
corresponding to the SLD integration (see solid curve) nearly 50 times! Note that
505
I.P.Omelyan, I.M.Mryglod, R.Folk
very small step sizes are impractical because then too much expensive force and
field evaluations should be done during the typical observation times. Moreover, the
ABM scheme appears to be about 2 times slower even if one iteration only is applied
within the corrector procedure.
Figure 1. The reduced total energy E∗/N (subsets (a)–(d)) and magnetization
S/N ((e)–(h)) per spin as functions of the length t/h of the simulations performed
on a Heisenberg ferrofluid within the decomposition SLD (solid curves), ABM
predictor-corrector (dashed curves in (a) and (b)), SLD4 (dashed curve in (e)),
and MTS (solid curves in (c) and (d)) algorithms at four fixed time steps: h∗ =
0.00125 ((a), (e)); 0.0025 ((b), (f)); 0.005 ((c), (g)); and 0.01 ((d), (h)). The initial
levels E∗(0)/N and S(0)/N are plotted by the thin horizontal line.
No systematic drift in E(t) and S(t) was observed within our decomposition
SLD algorithm at time steps up to h∗ = 0.01 over a length of t/h = 100 000.
The precision of the algorithm was measured in terms of the ratio ΓE = (〈(E(t)−
E(0))2〉/〈(U(t)−U(0))2〉)1/2 of total E and potential U energy fluctuations. Taking
into account that for our system (〈U(t) − U(0)〉2)1/2 ≈ 0.0335, we have obtained:
ΓE ≈ 0.12%, 0.28%, 0.98% and 7.7% for the time steps h∗ = 0.00125, 0.0025, 0.005
506
Spin dynamics simulations. . .
and 0.01, respectively. It is a common requirement, the ratio ΓE should not exceed
one or two per cent to properly reproduce features of microcanonical ensembles. As
we can see, time steps of h∗ < 0.01 satisfy this requirement and, thus, they can be
used for precise calculations. The case h∗ ∼ 0.01 can also be acceptable when the
precision is not so important, for example, in hybrid MC simulations [26], where
only the time reversibility is required to satisfy detailed balance.
Note that the decomposition and ABM methods lead to the total momentum
conservation within machine accuracy (∼ 10−12 in our program code). The reason
is that all velocities are updated simultaneously and the interparticle forces can be
evaluated exploiting the third Newton law. For similar reasons, the ABM integra-
tion maintains the total magnetization S (but it does not conserve individual spin
lengths). During the decomposed integration the magnetization will not be conserved
exactly. However, the fluctuations appear to be very small (see figure 1 (e)-(h)) and
consist 〈(S(t) − S(0))2〉1/2/N ≈ 10−7, 5 · 10−7, 2 · 10−6 and 10−5 at h∗ = 0.00125,
0.0025, 0.005 and 0.01, respectively.
An example on the total energy conservation obtained within the fourth-order
SLD4 integration (equations (7)–(10), (12) and (13)) at h∗ = 0.01 ≡ hmax is shown
in figure 1 (d). The energy fluctuations are decreased in this case approximately two
times, i.e., Γ
(4)
E (hmax) = ΓE(hmax)/2. This compensates to some extent the additional
computational efforts needed to evaluate high-order expressions (12). However, the
SLD4 algorithm allows to reduce the magnetization deviations to a negligible small
level which even for the greatest time step hmax = 0.01 does not exceed about 3·10−9.
Moreover, the time efficiency of the SLD4 algorithm will increase with decreasing of
the step size, because the energy fluctuations tend to zero in this case as Γ
(4)
E (h) ∼ h4
at h → 0, i.e., much faster than during the second-order SLD integration, where
ΓE(h) ∼ h2. The limiting value h̃ at which the integration will be performed with
the same efficiency both for SLD and SLD4 algorithms can be estimated using the
relation Γ
(4)
E (h̃) = ΓE(h̃/P ), where the fact that the five-stages SLD4 algorithm is
in factor P = 5 slower than has been used. Then in view of the explicit dependence
of ΓE and Γ
(4)
E on h and the above boundary condition Γ
(4)
E (hmax) = ΓE(hmax)/2
one obtains that h̃ = hmax
√
2/P ≈ 0.003. Therefore, the forth-order version SLD4
can be recommended if the total energy should be conserved with the precision
ΓE(h̃/P ) = ΓE(hmax)(h̃/(hmaxP ))2 = ΓE(hmax)2/P
4 ≈ 0.025% or better.
Finally, the function E(t) corresponding to the MTS integration (equations (9),
(11) and (15)) is plotted in figure 1 (c) and (d). The multistep parameter k (ap-
pearing in equation (15)) was chosen to be equal to 3. As can be seen, the MTS
algorithm conserves the total energy more accurately than the usual SLD integrator.
For instance at h∗ = hmax, the relative energy fluctuations are reduced in factor four,
i.e., Γmts
E (hmax) = ΓE(hmax)/4. To achieve such a reduction within the usual scheme,
it is necessary to decrease the step size h two times, and thus to increase the pro-
cessor time by 100% for covering the fixed interval of integration. Within the MTS
algorithm, the processor time needed to perform additional 2k evaluations of fast
varying forces consists of 45% only on the total computational time. Note that the
fast varying forces f
(f)
i were updated in our case in ∼ (Rc/R)3 = (2.5/2(1/6))3 ≈ 11
507
I.P.Omelyan, I.M.Mryglod, R.Folk
times faster than the long-ranged terms f
(s)
i . The efficiency of the MTS integration
will further increase (for larger systems) with increasing of the ratio R c/R.
5. Investigation of collective excitations
The dynamic behaviour of the Heisenberg fluid was investigated by measuring
the spectra F (k, ω) = 1
2π
∫∞
−∞
F (k, t)e−iωtdt of collective density-density Fnn(k, t),
spin-spin F L,T
ss (k, t), and spin-density F L
sn(k, t) time correlation functions,
Fnn(k, t) =
1
N
〈
N
∑
i,j
eik·(ri(0)−rj(t))
〉
,
F L,T
ss (k, t) =
1
N
〈
N
∑
i,j
s
L,T
i (0)· sL,Tj (t)eik·(ri(0)−rj(t))
〉
,
F L
sn(k, t) =
1
N
〈
N
∑
i,j
sLi (0)e
ik·(ri(0)−rj(t))
〉
≡ F L
ns(k, t) . (18)
Note that the superscripts (L) and (T) refer to the longitudinal and transverse
components of si with respect to the magnetization vector S rather than with respect
to k. The time correlation functions (18) were calculated using our SLD algorithm
at h∗ = 0.005, and the microcanonical averaging 〈 〉 was taken over a length of 100
000 steps. This length is sufficient for the highest frequency resolution of F (k, ω).
Ten independent runs with different initial equilibrium configurations were carried
out to reduce the statistical noise.
The wavevector-dependent correlation functions Fnn(k), F L
sn(k), FT
ss (k), and
F L
ss(k) obtained in the static limit (t = 0) are shown in figure 2. The density-density
structure factor Fnn(k) is related to a binary distribution function of magnetic par-
ticles averaged over their spin orientations. At intermediate and large wavevector
values of k this factor has a similar shape as that corresponding to pure liquid sys-
tems [27]. However, at small wavevectors we can observe deviations from the liquid
behaviour. For example, the function Fnn(k) has an additional minimum at kσ ∼ 1.5
(see subset (a)) that is caused by the existence of spin-spin interactions. The spin-
density structure factor F L
sn(k) varies on wavevector like Fnn(k), but it appears to
be shifted down at not very small k. This is so because in the infinite-wavevector
regime k → ∞ the function F L
sn(k) should tend to the magnetization S/N = 0.6536
of the system which cannot exceed unity, whereas limk→∞ Fnn(k) = 1 by definition.
At the same time, the transverse F T
ss (k) and longitudinal F L
ss(k) components of
the spin-spin static correlation function differ between themselves in a characteristic
way (see subset (b) of figure 2). While the transverse function tends to infinity at
k → 0 (this manifests itself the presence of Goldstone modes), the longitudinal
component in this regime accepts finite values related to the magnetic susceptibility
of the system. Note that the difference between F T
ss (k) and F L
ss(k) is due to the
fact that we are investigating a ferrophase. In the case of a paraphase when the
508
Spin dynamics simulations. . .
Figure 2. The static density-density Fnn(k) and spin-density F L
sn(k) (subset (a))
as well as transverse FT
ss (k) and longitudinal F L
ss(k) spin-spin (subset (b)) corre-
lation functions of a Heisenberg ferrofluid as depending on wavenumber.
magnetization is equal to zero (at the absence of an external field), these functions
will be identical.
The density-density Fnn(k, ω) and longitudinal spin-density F L
sn(k, ω) dynamic
structure factors are plotted in figures 3 and 4, respectively. We have used a nondi-
mensional presentation, F ∗(k, ω) = F (k, ω)/τtr, for all the spectral functions with
ω∗ = ωτtr and k∗ = k/kmin, where kmin = 2π/V 1/3 ≈ 0.53/σ is the smallest accessible
wavevector. The peak observed in Fnn(k, ω) at non-zero frequency ω = ωs(k) 6= 0 for
each wavevector value k 6= 0 should be associated with a propagator sound mode
well established for liquid systems [27,28]. This peak is especially visible in subset (a)
of figure 1 corresponding to the smallest wavevector k = kmin. At the same time,
a maximum of Fnn(k, ω) at ω = 0 represents the well-known diffusive heat mode.
Therefore, like pure liquids the density spectrum of the Heisenberg ferrofluid can be
expressed in terms of one heat and two complex-conjugate sound modes with the
dispersion ωs(k). The spin-density function F L
sn(k, ω) behaves similar to Fnn(k, ω).
Note that the transverse component FT
sn(k, ω) vanishes for arbitrary values of k and
ω because of 〈sTi 〉 = 0.
Consider, finally, the process of spin relaxation described by the wavevector- and
frequency-dependent transverse FT
ss (k, ω) and longitudinalF L
ss(k, ω) parts of the spin-
spin correlation function. These functions are shown in figures 5 and 6, respectively.
For the transverse function FT
ss (k, ω) we have identified one peak at ω = ω∗
T(k) for
each wavevector k considered. This peak is very sharp at small k (see subset (a)
of figure 5) and shifts to the right with increasing k (see subsets (b), (c) and (d)).
Like lattice systems, such a quasiparticle behaviour should be associated with the
existence of transverse spin waves in the Heisenberg ferrofluid. Up to three maxima
were observed (see figure 6) for the longitudinal function F L
ss(k, ω). In particular at
509
I.P.Omelyan, I.M.Mryglod, R.Folk
Figure 3. The density-density dynamic structure factor of a Heisenberg ferrofluid
as a function of reduced frequency ω∗ at four fixed values of reduced wavenumbers:
k∗ = 1 (subset (a)), k∗ = 2 (b), k∗ = 3 (c), and k∗ = 5 (d).
Figure 4. The spin-density dynamic structure factor of a Heisenberg ferrofluid as
a function of reduced frequency ω∗ at four fixed values of reduced wavenumbers:
k∗ = 1 (subset (a)), k∗ = 2 (b), k∗ = 3 (c), and k∗ = 5 (d).
510
Spin dynamics simulations. . .
Figure 5. The transverse spin-spin dynamic structure factor of a Heisenberg
ferrofluid as a function of reduced frequency ω∗ at four fixed values of reduced
wavenumbers: k∗ = 1 (subset (a)), k∗ = 2 (b), k∗ = 3 (c), and k∗ = 5 (d).
Figure 6. The longitudinal spin-spin dynamic structure factor of a Heisenberg
ferrofluid as a function of reduced frequency ω∗ at four fixed values of reduced
wavenumbers: k∗ = 1 (subset (a)), k∗ = 2 (b), k∗ = 3 (c), and k∗ = 5 (d).
511
I.P.Omelyan, I.M.Mryglod, R.Folk
k∗ = 1, the frequency positions for these maxima are: ω∗
L
(1) = 0, ω∗
L
(2) ≈ ω∗
T ≈ 1 and
ω∗
L
(3) ≈ ω∗
s ≈ 2.37. While the first maximum at ω = 0 corresponds to pure diffusive
relaxation processes, the position of the second maximum coincides with that of the
transverse spin wave peak, indicating the possibility for propagation of longitudinal
spin waves as well (which however are dumped more rapidly). The origin of the
third maximum can be explained by direct effect of the liquid subsystem on spin
relaxation, since its position coincides completely with the sound mode peak in the
density spectrum.
In general, the results obtained for dynamic structure factors are in good agree-
ment with predictions of [29]. The additional possibility of longitudinal spin wave
propagation in magnetic liquids at sound frequency can be considered as a new
effect which has yet to be quantitatively described in theory as well as observed
experimentally.
6. Conclusion
Starting from the Liouville operator formalism in tandem with Suzuki-Trotter-
like decompositions of exponential operators, we have devised and tested a class of
new algorithms for numerical integration of the equations of motion at the presence
of both liquid and spin subsystems. The greatest advantages of the new algorithms
over standard predictor-corrector schemes can be summarized as follows: (i) sym-
plecticity and time-reversibility, (ii) exact conservation of individual spin lengths,
(iii) explicitness (no iteration), (iv) much more accuracy in the total energy con-
servation. Moreover, their excellent stability makes it possible to apply much larger
time steps that may lead, therefore, to a substantial speedup of MD simulations on
magnetic liquids.
As an example of an interesting physical system we have studied a three-dimensi-
onal isotropic Heisenberg ferrofluid model. Using the algorithms proposed has al-
lowed one to determine quantitatively the density-density, spin-spin and spin-density
dynamic structure factors for the first time. As a result, we have shown that like
pure liquid systems the density spectrum can be described by one diffusive heat
and two propagative sound modes. At the same time similar to spin lattices, the
propagative transverse and longitudinal spin excitations are made of one primary
wave peak at each wavevector. In the spectrum of longitudinal spin excitations, one
secondary wave peak has been predicted additionally for the Heisenberg ferrofluid.
It has been identified that the frequency position of this peak is identical to that of
a sound mode. This result evidently follows from the effect of the liquid subsystem
on spin dynamics.
The presence of a secondary wave peak in magnetic liquids indicates the possi-
bility of longitudinal spin wave propagation at sound frequency. Such propagation
should be considered as a new effect which has yet to be directly tested experimen-
tally. The question of how to describe quantitatively this effect in theory will be a
subject of our further studies.
512
Spin dynamics simulations. . .
Acknowledgements
Part of this work was supported by the Fonds zur Förderung der wissenschaftli-
chen Forschung under Project No. P12422–TPH.
References
1. Amit D.J. Field Theory, the Renormalization Group, and Critical Phenomena. New
York, McGraw-Hill, 1978.
2. Handrich K., Kobe S. Amorphe Ferro- und Ferrimagnetika. Berlin, Akademika-Verlag,
1980.
3. Parisi G. Statistical Field Theory. Wokingham, Addison-Wesley, 1988.
4. Lomba E., Weis J.J., Almarza N.G., Bresme F., G. Stell. // Phys. Rev. E, 1994,
vol. 49, p. 5169.
5. Tavares J.M., Gama M.M.T., Teixeira P.I.C., Weis J.J., Nijmeijer M.J.P. // Phys.
Rev. E, 1995, vol. 52, p. 1915.
6. Nijmeijer M.J.P., Weis J.J. // Phys. Rev. Lett., 1995, vol. 75, p. 2887; Phys. Rev. E,
1996, vol. 53, p. 591.
7. Chen K., Landau D.P. // Phys. Rev. B, 1994, vol. 49, p. 3266.
8. Evertz H.G., Landau D.P. // Phys. Rev. B, 1996, vol. 54, p. 12302.
9. Bunker A., Chen K., Landau D.P. // Phys. Rev. B, 1996, vol. 54, p. 9259.
10. Costa B.V., Costa J.E.R., Landau D.P. // J. Appl. Phys., 1997, vol. 81, p. 5746.
11. Landau D.P., Krech M. // J. Phys. Cond. Mat., 1999, vol. 11, p. R1.
12. Ido Y., Tanahashi T. // J. Phys. Soc. Jap., 1991, vol. 60, p. 466.
13. Akhiezer I.A., Akhiezer I.T. // Sov. Phys. JETP, 1984, vol. 59, p. 68.
14. Akhiezer I.A., Akhiezer I.T. // Sov. Phys. Solid State., 1987, vol. 29, p. 48.
15. Mryglod I.M., Tokarchuk M.V., Folk R. // Physica A, 1995, vol. 220, p. 325.
16. Mryglod I.M., Folk R. // Physica A, 1996, vol. 234, p. 129.
17. Burden R.L., Faires J.D. Numerical Analysis. 5th ed. Boston, PWS Publishing, 1993.
18. Krech M., Bunker A., Landau D.P. // Comput. Phys. Commun., 1998, vol. 111, p. 1.
19. Swope W.C., Andersen H.C., Berens P.H., Wilson K.R. // J. Chem. Phys., 1982,
vol. 76, p. 637.
20. Frenkel D., Smit B. Understanding Molecular Simulation: from Algorithms to Appli-
cations. New York, Academic Press, 1996.
21. Allen M.P., Tildesley D.J. Computer Simulation of Liquids. Oxford, Clarendon, 1987.
22. Tuckerman M., Berne B.J., Martyna G.J. // J. Chem. Phys., 1992, vol. 97, p. 1990.
23. Suzuki M., Umeno K. – In: Computer Simulation Studies in Condensed Matter
Physics VI, edited by Landau D.P., Mon K.K., Schüttler H.B., Springer, Berlin, 1993.
24. Folk R., Mryglod I.M., Omelyan I.P. (unpublished).
25. Craknell R.F., Nicholson D., Parsonage N.G. // Mol. Phys., 1990, vol. 71, p. 931.
26. Duane S., Kennedy A.D., Pendleton B.J., Roweth D. // Phys. Lett. B, 1987, vol. 195,
p. 216.
27. Mryglod I.M., Omelyan I.P. // Mol. Phys., 1995, vol. 84, p. 235.
28. Mryglod I.M., Omelyan I.P. // Mol. Phys., 1997, vol. 91, p. 1005.
29. Mryglod I., Folk R., Dubyk S., Rudavskii Yu. // Physica A, 2000, vol. 277, p. 389.
513
I.P.Omelyan, I.M.Mryglod, R.Folk
Спiн-динамiчнi розрахунки колективних збуджень у
магнiтних piдинах
I.П.Омелян 1 , І.М.Мриглод 1,2 , Р.Фольк 2
1 Інститут фізики конденсованих систем НАН Укpаїни,
79011 Львів, вул. Свєнціцького, 1
2 Інститут теоретичної фізики, Університет м. Лінц,
Альтенбергштрассе 69, A-4040 Лінц, Австрія
Отримано 29 квiтня 2000 р.
Пpопонується новий пiдхiд до дослiдження динамiчних властивостей
спiнових piдин у комп’ютеpному експеpиментi. Пiдхiд базується на
фоpмалiзмi опеpатоpiв Лiувiлля для опису Гамiльтонової динамiки у
поєднаннi з методом фактоpизацiї експонентних пpопагатоpiв типу
Сузукi-Тpоттеpа. У pезультатi розвинуто сукупність нових алгоритмів
для числового iнтегpування piвнянь pуху за наявностi як тpасляцiй-
них, так i спiнових ступеней вiльностi, які зберігають фазовий об’єм
і є часо-зворотніми. Показано, що цi алгоpитми можуть викоpисто-
вуватися у комп’ютеpних розрахунках із набагато бiльшими часови-
ми кpоками порівняно з тими, якi пpитаманнi стандаpтним пpогнозо-
коpектованим схемам. Це нам дало змогу впеpше виконати прямі
розрахунки динамiчних стpуктуpних фактоpiв типу спiн-спiн, спiн-
густина i густина-густина для Гайзенбеpгiвської моделi феpофлюїду.
Встановлено, що подiбно до чистих piдин густинний спектp може бу-
ти описаний в термінах теплової та звукових мод. Водночас, подiб-
но до спінових гpаткових систем,у феромагнітній фазі у поздовж-
ньому та попеpечному динамічних структурних факторах спін-спін
як основний спостерігається спін-хвильовий пік. Окрім того, у поз-
довжньому динамічному структурному факторі спін-спін виявлено
додатковий пік, поява якого пеpедбачалася нами раніше [Mryglod I.,
Folk R. et al., Physica A277 (2000) 389]. Частотна позиція цього піку
повністю збігається з частотою звукових збуджень, що відображає
вплив рідинної підсистеми на спінову динаміку. Можливiсть поши-
pення поздовжнiх спiнових хвиль на звуковій частоті слід розглядати
як новий ефект, що потребує експериментальної перевірки.
Ключові слова: Гайзенбеpгiвська piдина, спiнова pелаксацiя,
колективнi моди, симплектичне iнтегpування, фактоpизацiї
Сузукi-Тpоттеpа
PACS: 02.60.Cb, 75.40.Gb, 75.50.Mm, 76.50.+g, 95.75.Pq
514
|