PHYSICAL REVIEW B 70, 064401 (2004)

Dzyaloshinskii-Moriya interactions and adiabatic magnetization dynamics in molecular magnets
H. De Raedt,1,* S. Miyashita,2,† K. Michielsen,1,‡ and M. Machida2,§

1Department of Applied Physics, Materials Science Centre, University of Groningen, Nijenborgh 4,

NL-9747 AG Groningen, The Netherlands
2Department of Applied Physics, Graduate School of Engineering, University of Tokyo, Bunkyo-ku, Tokyo 113-8656, Japan(Received 23 December 2003; revised manuscript received 27 April 2004; published 3 August 2004)
A microscopic model of the molecular magnet V15 is used to study mechanisms for the adiabatic change ofthe magnetization in time-dependent magnetic fields. The effects of the Dzyaloshinskii-Moriya interaction, themost plausible source for the energy-level repulsions that lead to adiabatic changes of the magnetization, arestudied in detail. We find that the energy-level repulsions that result from this interaction exhibit a strong
dependence on the direction of the applied field. We also discuss the role of magnetic anisotropy in themolecule Mn12-acetate.
DOI: 10.1103/PhysRevB.70.064401

PACS number(s): 75.10.Jm, 75.50.Xx, 75.45.!j, 75.50.Ee

I. INTRODUCTION

Recently, magnetic molecules such as Mn12 or V15 haveattracted a lot of interest. These nanomagnets are often used
to study explicit real-time quantum dynamics, e.g., tunnelingof the magnetization and quantum (de)coherence.1–22 As aresult of the very weak intermolecular interactions, experiments can directly probe the magnetization dynamics of theindividual molecules. In particular, the adiabatic change ofthe magnetization at low temperature is governed by the discrete energy-level structure.23–27
The adiabatic change of the magnetization requires someinteractions that yield energy-level repulsions, i.e., interactions that do not commute with the magnetization. The
Dzyaloshinskii-Moriya (DM) interaction is the most likely
candidate for such an interaction.11,28,29 In the case of anisotropic high-spin molecules such as Mn12 and Fe8, simplified
anisotropic single-spin models for a specific spin multiplet
can approximately reproduce the gaps of the level repulsions.
However, the case of V15 is more complicated because V15
has half-odd-integer spin and the time-reversal symmetry enforces at least a twofold degeneracy of the energy levels atzero field.
The possibility that the DM interaction might be the main
mechanism for the adiabatic change for the magnetization in
V15 has been explored in earlier work29,30 but the dependenceof the energy-level scheme on the direction of the magnetic
field was not considered. In this paper we further elaborateon this issue and demonstrate that the effect of the DM interaction on the magnetization dynamics strongly depends onthe direction of the applied magnetic field. This directional
dependence has not been observed in experiments on V15.
Therefore the DM interaction generally does not explain why
the magnetization changes as the magnetic field is swept.
It has been pointed out that the DM interaction is accompanied by a higher-order correction term that restores theSU(2) symmetry.31–34 In this paper, we focus on the effectsof the DM interaction and leave the inclusion of the higherorder term for future study. As another source of level repulsion, we might consider the hyperfine interaction with thenuclear spin. The effects of the hyperfine interaction havealready been discussed in Ref. 52.
1098-0121/2004/70(6)/064401(12)/$22.50

As the DM interaction has a vector character and is anisotropic, the dynamics of the magnetization is expected todepend on the direction of the magnetic field. First we study
the characteristic properties of the DM interaction for a simplified model of V15, namely, three spins on a triangle. Then
we confirm the properties found in the three-spin model byfull diagonalization of the 15-spin model of V15. In order tobridge the energy scales involved (e.g., from 800 K, a typicalenergy scale for the interaction between individual magnetic
ions, to about 10−2 K, a typical energy scale for energy-level
splittings), the calculation of the energy levels of the manyspin Hamiltonian has to be very accurate. We have tested
various standard algorithms to compute the low-lying states.
For systems that are too large to be solved by full exact
diagonalization (such as the 15-spin V15 model), we use theLanczos method with full orthogonalization (LFO), a Chebyshev polynomial projector (CPP) method, and a power
method with additional subspace diagonalization. These algorithms can solve the rather large eigenvalue problems with
sufficient accuracy. The consistency of the data obtained bydifferent methods gives extra confidence in the numericalresults.
The magnetic properties of molecules such as Mn12 areoften studied by considering a simplified model for the magnetic energy levels for a specific spin multiplet, e.g., S = 10.
However, for these and other, similar, magnetic moleculesthat consist of several magnetic moments [in the case of
Mn12, eight Mn3+ !S = 2" and four Mn4+ !S = 3 / 2"], the reduction of the many-body Hamiltonian to an effective Hamiltonian for a specific spin multiplet is nontrivial. Magnetic
anisotropy, a result of the geometrical arrangement of themagnetic ions within a molecule of low symmetry, mixes
states of different total spin and enforces a treatment of thefull Hilbert space of the system. For Mn12, the dominant
contribution to the magnetic mixing due to spin-orbit interactions is also given by the DM interaction.35,36 In principle,
this type of interaction can change energy-level crossings
into energy-level repulsions. The presence of the latter isessential to explain the adiabatic changes of the magnetization at the resonant fields.23–27 Thus, a minimal magnetic
model Hamiltonian should contain (strong) Heisenberg inter-

70 064401-1©2004 The American Physical Society

PHYSICAL REVIEW B 70, 064401 (2004)

De RAEDT et al.

actions, DM interactions, and a coupling to the applied magnetic field.11,28–30,37–42 Experiments on Mn12 suggest that theenergy gaps related to the transition from a state with magnetization M # −10 to a state with magnetization M " 4 areof the order of 10−9 K.43 Such gaps are too small to detectwith standard precision (13–14 digits) calculations, andtherefore in this paper we present only the global energylevel diagram obtained from microscopic model calculations.
The paper is organized as follows. In Sec. II we analyze areduced three-spin model with C3 symmetry for the V15 molecule. Results for the energy-level schemes for a 15-spin
model of the V15 are presented in Sec. III. In Sec. IV we
discuss the effects of anisotropic terms and the cases withless symmetry. In Sec. V, we report results for Mn12. In Sec.
VI we give our conclusions. In Appendixes A and B, webriefly discuss some analytical solutions and the numericalalgorithms that we use to compute the energy levels, respectively.
II. TRIANGLE MODEL WITH C3 SYMMETRY

The V15 molecule has C3 symmetry but is not invariant
for mirror reflection about the triangle plane. Therefore wetake as a simplified model for the V15 molecule, a system ofthree spins on a triangle with C3 symmetry only. We choose
the z axis to lie along the axis of C3 symmetry. The Hamiltonian is given by
Ji,i+1Si · Si+1 + $ Di,j · 'Si # S j( − h · ) $ Si* ,
$
i=1
%i,j&
i
3

H=−

FIG. 1. (Color online) Energy levels of Hamiltonian (1) for J= −2.5 K, Dx = Dy = Dz = 0.25 K. At each value of h the slope of eachlevel gives the corresponding value of the projection of the totalmagnetization on the magnetic field axis. Top: Applied magnetic
field h parallel to the z axis. Bottom: Applied magnetic field h along
the x axis.

!1"
where J1,2 = J2,3 = J3,1 + J denotes the exchange interaction
and h represents the applied magnetic field. In general we
can choose any direction of the DM vector Di,j unless there
is some additional symmetry. In the present case, because of
the C3 symmetry, the z component of the DM vectors mustall be equal, i.e.,
z
z
z
D1,2
= D2,3= D3,1+ Dz .

!2"

If the system has mirror reflection symmetry about the triangle plane then Dz = 0. The x and y components of the DMy, must obey the relation
vector, Dxi,j and Di,j

, .
) *,
.) *
) *
x
D2,3
y
D2,3

x
D3,1
y
D3,1

1
2

−

=

−

−

=

-3
2

-3
2

1
−
2

-3
1
−
2
2

-3
2

x
y
where Dx + D1,2
and Dy + D1,2
.

1
−
2

) *
Dx

Dy

the applied field is parallel to the z axis and the conditions(2) and (3) hold. Then the Hamiltonian (1) is block diagonal,
the matrix containing four blocks of 2 # 2 matrices. The(eight-dimensional) Hilbert space separates into four twodimensional spaces:

,

Dx
,
Dy

!3"

In Appendix A we give the analytic expressions for the
eigenvalues and eigenvectors of model (1) in the case that/03/2,3/2&,0a&1,/03/2,1/2&,0ā&1,/03/2,1/2&,0b̄&1,/03/2,3/2&,0b&1,!4"

where the state 0S , M& denotes an eigenstate of the Heisenberg model with total spin S and magnetization M. The expressions for the orthonormal states 0a& and 0b& are given in
Appendix A. The state 0ā& !0b̄&" denotes the state 0a& !0b&"with all spins reversed. From Eq. (4) it follows that there isno mixing among the four levels with S = ± 1 / 2, and therefore they simply cross each other. From the analytical solution it is also easy to see that there are no energy-level repulsions at h = 0. Furthermore, one can show analytically that
it is impossible to change the magnetization from −3 / 2 to
+ 3 / 2 by (adiabatically) reversing the external field along thez axis.28
In the top panel of Fig. 1 we show the energy levels as afunction of the strength h of the applied magnetic field forthe case Dx = Dy = Dz. The field is aligned along the z axis. It
is evident that there is no level mixing at h = 0. It is important
to note that the small energy difference between the first

064401-2

PHYSICAL REVIEW B 70, 064401 (2004)

DZYALOSHINSKII-MORIYA INTERACTIONS AND…

FIG. 2. (Color online) Top: Energy levels of Hamiltonian (1) forJ = −2.5 K, Dx = Dy = Dz = 0.25 K, and the applied magnetic field h= h!1 , 0 , 1" / -2 tilted by 45° with respect to the z axis. Bottom:
Detailed view of the h dependence of the four lowest energy levels.
At each value of h the slope of each level gives the corresponding
value of the projection of the total magnetization on the magnetic
field axis.

excited state and the degenerate ground states at h = 0 is not
the minimum energy difference of the avoided level crossing(near h = 0). Therefore, at h = 0 there is no adiabatic change ofthe magnetization. On the other hand, if we apply the magnetization along the x direction (i.e., perpendicular to thesymmetry axis of the triangle), then two of the states become
degenerate. We find two nearly degenerate avoided levelcrossings, as shown in the bottom panel of Fig. 1. As anintermediate case, in Fig. 2 we show the energy-level diagram for the case that the magnetic field is tilted by 45° withrespect to the z axis 'h = h!1 , 0 , 1" / -2(. Then, there is a
simple crossing at zero field and an avoided level crossing
between the levels of M = ± 1 / 2 at a finite magnetic field.
Indeed, a closer look at the level diagram (see the bottom
panel in Fig. 2) reveals that the minimum energy difference
between the two pairs of levels does not occur at zero field
but at h # 0.35 T. This implies that the Landau-ZenerStückelberg transition from the 01 / 2 , −1 / 2& to the 01 / 2 , 1 / 2&level does not take place at h = 0 but at h # 0.35 T. In conclusion, the position and the energy splitting of the avoided
level crossing that is responsible for the adiabatic change
depend on the direction of the field.
The numerical results discussed above have been obtainedfor Dx = Dy = Dz. In Ref. 22 the DM vector is taken parallel tothe y axis at all the bonds and the field is applied along the zaxis. For the present model, this case corresponds to the casewith only Dz !Dx = Dy = 0" and a field applied in the x direc-

FIG. 3. Schematic diagram of the magnetic interactions inmodel (1) of the V15 molecule.

tion. In this case, the gap opens symmetrically with respect
to the applied field.22 However, Ref. 22 did not address thedependence on the direction of the magnetic field. In the next
section we repeat the analysis of this dependence for a 15spin model of the V15 molecule.
III. 15-SPIN MODEL FOR THE VANADIUM COMPLEX
V15
A. Spin interactions in V15

In Fig. 3, we show the schematic diagram of the dominantmagnetic (Heisenberg) interactions in the 15-spin model ofIV
As6O42!H2O"( · 8H2O). The magthe V15 molecule (K6'V15
netic structure consists of two hexagons with six S = 1 / 2
spins each, enclosing a triangle with three S = 1 / 2 spins. All
dominant Heisenberg interactions are antiferromagnetic. The
dimension of the Hilbert space of this model is 215 = 32 768.
The minimal Hamiltonian is given by expression (1) with 15
instead of three spins.22,28,29,40 The Heisenberg interactions
Ji,j in Eq. (1) between the vanadium atoms are defined according to Fig. 3. For simplicity, we assume that Di,j = 0 except for bonds for which the Heisenberg exchange constantis the largest (i.e., equal to J).29,40 Rotations of 2$ / 3 and
4$ / 3 around the axis perpendicular to and passing through064401-3

PHYSICAL REVIEW B 70, 064401 (2004)

De RAEDT et al.

the center of the hexagons leave the V15 complex invariant.
This enforces the constraints (2) and (3) on the values ofDi,j.29,30
B. Energy-level diagrams

In this section, we will study the energy eigenvalues forvarious sets of parameters of the Hamiltonian, i.e., Jij and/orDij, which have so far been proposed for the 15-spin modelof V15.4,18,29,40
First let us study the case of the model parameters given
in Ref. 40,
J = − 800 K,

J1 = J! = − 54.4 K,

J2 = J" = − 160 K,!5"J3 = J4 = J5 = J6 = 0.

In the absence of DM interactions, our method reproduces
the energy gap between the ground state and the first excitedstate at h = 0 to be a value of 4.124 78 K in agreement withthe result of Ref. 40. This value is in reasonable agreementwith the experimental value of approximately 3.7 K.43 This
energy gap causes a transition between the states 01 / 2 , 1 / 2&
and 03 / 2 , 3 / 2& which takes place at h # 2.8 T (because 1 T
corresponds to 1.343 K), in good agreement with the experimental value 2.8 T.
Following Ref. 30, we now take for the DM interaction
parametersxyzx
D10,11= −14.641 K, D10,11= 54.641 K, D10,11= 40 K, D12,13yzx= −40 K, D12,13 = −40 K, D12,13 = 40 K, and D14,15
yz= 54.641 K, D14,15 = −14.641 K, D14,15 = 40 K. For this choiceof model parameters, the eight lowest energies of the V15
model (1) for two values of the applied magnetic field (h= 0 and h = 4 T) along the z axis can be found in Table I.
From Table I we see that for zero field the DM interaction
splits the doubly degenerate doublets of S = 1 / 2 states into
two doublets of S = 1 / 2 states. The difference in energy between the doubly degenerate, first excited states and the two
fold degenerate ground states has a value of 0.0085 K. Thisvalue is much smaller than the experimental estimate of
0.05 K,22 but of the same order of magnitude as the value
cited in Ref. 29. The next four higher levels are S = 3 / 2
states. The energy-level splitting between the S = 3 / 2 and S= 1 / 2 states is 4.1 K.
Next, we study the parameter set of Ref. 29:

J = − 800 K,
J2 = J" = − 350 K,

which is approximately 5% of the largest Heisenberg coupling. Using the rotational symmetry of the hexagon it folxy= 14.641 K, D3,4
= −54.641 K,lows from Eq. (3) that D3,4zxyz
D3,4 = 40 K and D5,6 = −54.641 K, D5,6 = 14.641 K, D5,6
= 40 K. If the two hexagons are not equivalent we cannot
reduce the number of free parameters by using a symmetryargument (see below). However, for simplicity, we may assume that the !x , y" positions of the spins on the lower hexagons differ from those on the upper hexagon by a rotation
about $ / 3. This yields for the remaining model parameters
TABLE I. The eight lowest eigenvalues Ei and total spin Si of
the corresponding eigenstates of the V15 model (1) with model parameters taken from Ref. 40 for two values of the external appliedfield h parallel to the z axis. The values of the DM vectors are givenin the text. The distance between Ei and the exact eigenvalue closest
to Ei is %i = %&i 0 !H − Ei"2 0 &i&1/2 " 10−9 for i = 1 , … , 7.

J3 = J4 = J5 = J6 = 0.

!7"

In the absence of DM interactions, the energy gap between
the fourfold degenerate ground state and the first excitedstate is 3.61 K, in full agreement with the result of Ref. 29.
Note that this value of the gap is fairly close to the experimental value of 3.7 K.43 Taking for the DM interactionsx
x
D1,2= D14,15= 25 K,

!6"

x
y
z
D1,2
= D1,2= D1,2= 40 K,J1 = J! = − 225 K,

x
x
x
x
D3,4
= D5,6= D10,11= D12,13= − 12.5 K,yyyy= − D5,6= − D10,11= D12,13= − 21.5 K,
D3,4

!8"

our calculation for the splitting between the two doubly degenerate S = 1 / 2 levels yields 0.0037 K, about a factor of 2
larger than the value cited in Ref. 29. For the energy splitting
between the S = 1 / 2 and S = 3 / 2 levels we obtain 3.616 K
instead of the value 3.618 K given in Ref. 29. These differences seem to suggest that a perturbation approach29 for the
DM interaction may not be sufficiently accurate for quantitative purposes.44
Finally, we study a parameter set obtained from a firstprinciples calculation18
J = − 809 K, J! = − 120 K,
J" = 120 K, J1 = − 30 K, J2 = − 122 K, J3 = − 3 K,
!9"i

Ei !h = 0"

Si !h = 0"

Ei !h = 4 T"

Si !h = 4 T"

J4 = − 11 K, J5 = − 3 K, J6 = − 2 K

0
1
2
3
4
5
6
7

−3679.53623744−3679.53623744−3679.52777009−3679.52777009−3675.42943612−3675.42943612−3675.42325141−3675.42325141

0.51
0.51
0.51
0.51
1.50
1.50
1.50
1.50

−3683.51181131−3682.21997451−3682.18488706−3678.11784886−3676.84225573−3676.83951808−3672.74011178−3667.37940477

1.50
0.51
0.53
1.50
0.52
0.51
1.50
1.50

(see Table I in Ref. 18). This set yields an energy gap of
4.915 K, as given in Ref. 18. Adding a DM interaction withthe parameters (8) does not significantly change the energy
gap between the singlet and triplet states.
In Fig. 4 we show the results for the set (7) and the DM
interaction parameters of (6). For the energy gap at zero
field, we find 3.7 K (instead of 4.1 K for the J’s of Ref. 40),in good agreement with the experimental estimate of 3.7 K.43
The results for the zero-field energy gap suggest that there
are many different sets of model parameters that approxi064401-4

PHYSICAL REVIEW B 70, 064401 (2004)

DZYALOSHINSKII-MORIYA INTERACTIONS AND…

FIG. 4. (Color online) Top: The eight lowest energy levels of
V15 model (1) with model parameters taken from Ref. 29 as afunction of the applied magnetic field h parallel to the z axis. The
values of the DM vector are given in the text. Bottom: Detailed
view of the four lowest energy levels at h # 0.

mately reproduce the experimental gap between the singlet
and triplet states. However, as we show below, the energy
gap at zero field does not necessarily correspond to the gapof a level repulsion that is required for an adiabatic change ofthe magnetization.
We now study the energy-level diagram as a function of
the direction of the magnetic field. As shown in the bottom
panel of Fig. 4, when the field is parallel to the symmetry
axis the energy levels of M = ± 1 / 2 simply cross, just as in
the case of the top panel of Fig. 1. In Fig. 5, we present
results for the cases where the angle between the appliedfield and the z axis is 45° and 90°, respectively. Clearly, we
find the same type of dependence of the energy levels on theangle of the field as in the case of the three-spin model (seebottom panel of Fig. 2 and bottom panel of Fig. 1, respectively). Exactly the same qualitative features are obtained forthe other sets of parameters (7) and (9) (results not shown).
Up to now, we used DM vectors that satisfy the rotational
symmetry of a hexagon, and we took the same DM vectors
for the other hexagon for simplicity.29 However, if there issome symmetry that connects the upper and lower hexagons,
we have a relation between the DM vectors on both hexagons. In concert with the relations between the exchange
couplings, let us assume that the upper and lower hexagons
are related to each other by a 180° rotation around a vector
that passes through V atom number 7 (see Fig. 3) and the
middle of the line connecting the two other V atoms of thetriangle. This symmetry operation X transforms the sites

FIG. 5. (Color online) Same as the right panel in Fig. 4 exceptthat the applied magnetic field h = h!1 , 0 , 1" / -2 is tilted by 45° withrespect to the z axis (top) and h is along the x axis (bottom).

/!1 , 2" , !3 , 4" , !5 , 6"1 into /!14, 15" , !10, 11" , !12, 13"1. If
we place the y axis along the line through V atom number 7
and through the middle of the line connecting V atoms 8 and
9, and if we take D1,2 = !Dx , Dy , Dz" as the reference DM
vector, the other DM vectors are given by
D3,4 = R2D1,2,

D5,6 = RD1,2 ,D14,15 = XD1,2,

D10,11 = RXD1,2 ,D12,13 = R2XD1,2 .

!10"Here R denotes a rotation of the hexagons around 2$ / 3 inthe plane of the hexagons. The energy-level diagrams for theset (5) and (6), with the additional constraint imposed by theabove symmetry, are qualitatively similar to those obtainedin Figs. 4 and 5 (results not shown).
Summarizing, as in the case of the three-spin model, all
our results for the 15-spin V15 model clearly demonstratethat the mixing of levels strongly depends on the direction of
the magnetic field. It seems therefore that this dependence isa generic feature of the DM interaction.

064401-5

PHYSICAL REVIEW B 70, 064401 (2004)

De RAEDT et al.

energy-level diagrams are qualitatively similar to those of theundistorted triangle (see Figs. 1 and 2), that is, crossings at
h = 0 and level repulsions close to h = 0 but not at h = 0. However, in the case of an undistorted triangle with anisotropic,
antiferromagnetic exchange interactions (i.e., different forthe x , y, and z components of the spins) and Dx = Dy = Dz
= 0.25 K, the energy-level diagrams are qualitatively similar
to the ones of the undistorted triangle (results not shown).
V. MANGANESE COMPLEX: Mn12

FIG. 6. (Color online) Energy-level diagram of the Hamiltonian(1) of a distorted triangle with model parameters J1,2 = −2.5 K,
J2,3 = −2.0 K, J3,1 = −3.0 K, and Dx = Dy = Dz = 25 K. The appliedmagnetic field h = h!1 , 0 , 1" / -2 is tilted by 45° with respect to the z
axis.
IV. EFFECTS OF LOWER SYMMETRY

We now discuss the effects of distortion of the triangle
and anisotropic exchange interactions in the triangle modelof V15 on the energy-level diagram for h # 0. When the triangle is distorted (J1,2 ! J2,3 ! J3,1 ! J1,2) the degeneracy of
the two doublets at h = 0 is lifted, even if Dx = Dy = Dz = 0. In
Figs. 6 and 7 we show data for J1,2 = −2.5 K, J2,3 = −2.0 K,
and J3,1 = −3.0 K and Dx = Dy = Dz = 0.25 K for the appliedmagnetic field h = h!1 , 0 , 1" / -2 tilted by 45° with respect tothe z axis, and for field directions parallel to the z axis andalong the x axis. Unless the field is parallel to the z axis, theThe four inner Mn4+ ions in the Mn12 molecule[Mn12!CH3COO"16!H2O"4O12 · 2CH3COOH· 4H2O],
have
spin S = 3 / 2. The other eight Mn3+ ions have spin S = 2. The
dimension of the Hilbert space of this system is 44 # 58
= 108. If the total magnetization is a conserved quantity, it
can be used to block-diagonalize the Hamiltonian, allowing a
numerical study of models of this size.41,45 However, tostudy the adiabatic change of magnetization, we have to take
into account all the states, and the dimension of the matrix
becomes prohibitively large. Thus we need to simplify themodel in order to reduce the dimension. A drastic reductionof the number of spin states can be achieved by assumingthat the strong antiferromagnetic Heisenberg interaction J!
between an inner ion and its outer neighbor allows the replacement of the magnetic moment of an inner ion by aneffective S = 1 / 2 moment. The schematic diagram of thissimplified (but still complicated) model37 is shown in the top
panel of Fig. 8. The dimension of the Hilbert space of thismodel is 24 # 54 = 104. In the following we study this simplified model.
The Hamiltonian for the magnetic interactions of the simplified Mn12 model can be written as37

)$ *
4

H=−J
+

Si

i=1

$

!i,j"!B

FIG. 7. (Color online) Same as Fig. 6 except that the appliedmagnetic field h is parallel to the z axis (top) and along the x axis(bottom).

4

2

− J!

z
!Si+4
"2
$ Si · S j − Kz $i=1!i,j"!B8

Di,j · 'Si # S j( −h · Si ,
$
i=1

!11"

where the index 1 ' i ' 4 !5 ' i ' 8" refers to S = 1 / 2 !S = 2"
spins and B denotes the set of pairs B= /!1 , 5" , !1 , 8" , !2 , 5" , !2 , 6" , !3 , 6" , !3 , 7" , !4 , 6" , !4 , 8"1. The
first two terms describe the isotropic Heisenberg exchange
between the spins. The third term (Kz) describes the singleion easy-axis anisotropy of S = 2 spins. The fourth term represents the antisymmetric DM interaction in Mn12. The vector Di,j determines the DM interaction between the ith S= 1 / 2 spin and the jth S = 2 spin. The last term describes theinteraction of the spins with the external magnetic field h.
Note that the factor g(B is absorbed in our definition of h.
Model (11) reproduces experimental data, such as the splitting of the neutron scattering peaks, the results of electron
paramagnetic resonance measurements, and the temperaturedependence of the magnetic susceptibility.37
The first three terms in the Hamiltonian (11) conserve the8
Szi . The DM interacz component of the total spin M z = $i=1
tion, on the other hand, mixes states with different total
spins. Hence, the DM interaction can change level crossings

064401-6

PHYSICAL REVIEW B 70, 064401 (2004)

DZYALOSHINSKII-MORIYA INTERACTIONS AND…

its). Clearly there are states with total spin eight, nine, andten within these 21 lowest eigenstates. Although the totalmagnetization is not a good quantum number, we can label
the various eigenstates by their (calculated) magnetization.
The S = 10 single-spin model for Mn12
H = − D!Sz"2 − hSz ,!12"

where D denotes the uniaxial anisotropy, is often used as astarting point to interpret experimental results.6,7,12–14,40 The
energy levels of this model exhibit crossings at the resonant
fields h = ± Dn for n = −10, … , 10, in qualitative agreementwith our numerical results (shown in the bottom panel of Fig.
8) for the microscopic model (11). A fit of the first eight levelcrossings of model (12) to the data of Fig. 8 yields D# 0.74 K, in good agreement with experiments.6,7 The
Hamiltonian of the single-spin model (12) commutes withthe magnetization Sz and therefore its energy diagram displays only level crossings, no level repulsions. Adding anisotropy of the form C!S+4 + S−4 " to model (12) changes the estimated value of D and leads to level repulsions when themagnetization changes by 4.37,46–49
It is also of interest to compare the level splitting at h
= 0 obtained by lowest-order degenerate perturbation theoryof model (12) with fourth-order anisotropy of the form S+4+ S−4 (Refs. 47, 50, and 51)%El+1,l = 32D
FIG. 8. Top: Schematic diagram of the magnetic interactions of
the simplified model (11) of the Mn12 molecule (Ref. 37). Blackcircles, S = 1 / 2; open circles, S = 2. Also shown are the DM vectors(for i " j and Di,j = −D j,i). Bottom: (Color online) The 21 lowest
energy levels of the Mn12 model (11) as a function of the appliedmagnetic field h. Solid lines, eigenstates with S # 10; dashed lines,
eigenstates with S # 9. The applied magnetic field h is parallel tothe z axis.

into level repulsions. Therefore, the presence of the DM interaction may be sufficient to explain the experimentally observed adiabatic changes of the magnetization.
The fourfold rotational-reflection symmetry (S4) of the
Mn12 molecule imposes some relations between the DM vectors. There are only three independent DM parameters:37xyz
Dx + D1,8, Dy + D1,8
, and Dz + D1,8
, as indicated in Fig. 8.
The parameters of the model (11) have been estimated bycomparing experimental and theoretical data. In this paper
we will use the parameter set B from Refs. 37 and 42: J= 23.8 K, J! = 79.2 K, Kz = 5.72 K, Dx = 22 K, Dy = 0, and Dz= 10 K. Although the amount of available data is not sufficient to fix all these parameters accurately, we expect that thegeneral trends in the energy-level diagram will not change
drastically if these parameters change relatively little.
In Table II we present results for the energy and total spinof the 21 lowest states for h = 0 and h = 5 T (in these calculations h is parallel to the z axis). The numerical results
obtained by full exact diagonalization (LAPACK), the Lanczos
method with full orthogonalization (see Appendix B), andthe Chebyshev-polynomial-based projector method (see Appendix B) are the same to working precision (about 13 dig-) *1
!S + m"!
C m/2
,
16D
'!m/2 − 1" ! (2 !S − m"!

!13"

for m even (%El+1,l = 0 for m odd) with the result of the
microscopic model calculation based on model (11). In Eq.
(13) l denotes the perturbed eigenstates in increasing order ofenergy and m is the magnetic quantum number of the unperturbed states.47 Using the values D # 0.69 K, C / D # 5.7# 10−5 obtained by fitting the single-spin model to experimental data,47 the energy gap for m = 6 is given by%E13,12 # 0.000 22 K.

!14"Taking into account that (because of the presence of S # 9
states) the corresponding S # 10 levels of model (11) are the14th and 15th lowest energy levels, Table II shows that for
h = 0, %E13,12 = 0.000 42 K. In view of the uncertainties on
the estimates of the various model parameters, the differenceof only a factor of 2 is remarkably small. From this comparison, we may conclude that the DM interaction leads to energy gaps that are of the same order of magnitude as the gaps
due to the fourth-order terms S+4 + S−4 in the single-spin model.
In the bottom panel of Fig. 8 we show the results for thelowest 21 energy levels of the Mn12 model as a function of
the applied magnetic field as obtained by LFO. The appliedmagnetic field is parallel to the z axis. In Fig. 8 solid(dashed) lines represent eigenstates with S # 10 (9) (within
an error of about 10%). Also eigenstates with S # 8 appear
for h ) 4 but these are not shown for clarity. For the Mn12
model, the DM induced energy splittings between the S# 10, M # −10 state and other states are less than 10−6 K.
Adding an extra transverse field by tilting the h field by 5°with respect to the z axis does not change this conclusion(results not shown).

064401-7

PHYSICAL REVIEW B 70, 064401 (2004)

De RAEDT et al.

TABLE II. The 21 lowest eigenvalues Ei and total spin Si of the corresponding eigenstates of the Mn12
model (11) for two values of the external applied field h along the z axis. The distance between Ei and theexact eigenvalue closest to Ei is %i = %&i 0 !H − Ei"2 0 &i&1/2 " 10−10 for i = 1 , … , 7. Note that for h = 0 the levels
12, 13, 18, 19, and 20 belong to the S # 9 subspace and not to the S # 10 subspace.
i

Ei !h = 0"

Si !h = 0"

Ei !h = 5 T"

Si !h = 5 T"

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

−812.771882673675−812.771882673460−798.326618260922−798.326618261122−785.677644659194−785.677644658983−774.774953284432−774.774953284294−765.549187817101−765.549187333902−757.919915510036−757.919915509970−757.673722613912−757.673722613981−751.806498916496−751.806072514140−747.135398548595−747.135398548602−746.357522623039−746.357522623082−745.778951523327

9.72
9.72
9.72
9.72
9.70
9.70
9.68
9.68
9.65
9.65
9.61
9.61
8.77
8.77
9.57
9.57
9.54
9.54
8.77
8.77
8.70

−878.203468556749−857.042137859145−837.727273846391−820.218084451590−816.530388063144−804.449056631056−804.242954124067−798.519706376890−790.336285910398−788.717912654407−782.037906282298−777.789121186874−774.222534072884−773.614730624955−766.993852023410−766.720405060281−760.807893152451−760.482768227423−757.060762637193−754.700878489864−753.310517350023

9.77
9.78
9.76
9.73
8.82
9.69
8.77
8.82
9.65
8.78
8.82
9.60
8.79
8.80
8.80
9.55
8.79
8.12
9.51
8.59
8.78

Finally, we added to model (11) the next-to-lowest orderrelativistic correction to the local anisotropy that is compatible with the symmetry of the square37
H1 = K1'!Sx1"2 + !S2y "2 + !Sx3"2 + !S4y "2(.

!15"

Although we took a perhaps unrealistically large value of
K1 !K1 = Kz / 2", we were unable to detect energy-level repulsions up to the M # −10, M # 3 transition (results notshown). On the other hand, in experiments,6,7,12 adiabaticchanges of the magnetization have been observed at h
# 3.4 T !M z # −10→ M z # 4"
and h # 3.9 T !M z # −10→ M z # 3" and the magnitude of the energy splittings is of
the order of 10 nK.43 The precision of the present calculations is about 10−6 K. Thus, it is consistent that within the(very high) resolution in the h field and the 13-digit precisionof the calculation, no information about the gap could be
extracted. The algorithms developed for the work presented
in this paper can be used for 33-digit calculations without
modification.
It has already been demonstrated that the many-spin
model (11) can reproduce neutron scattering data, highfrequency electron spin resonance, and the temperature dependence of the magnetic susceptibility of the Mn12
system.37 Our results for the energy-level diagram show features (for example, the presence of S = 9 levels, see Fig. 8)that are not accounted for by a single-spin S = 10 model for

Mn12. In Sec. III we discussed the properties of the level
repulsions due to the DM interactions in the V15 system. For
the Mn12 system we expect to find similar behavior. However, we find that the energy differences at the crossing
points are smaller than the 10−6 K accuracy of our numericalcalculations. Therefore in this section we studied the mainfeatures of the energy-level diagram of Mn12. In principle, it
is possible to leave the S # 10 multiplet by sweeping themagnetic field but for Mn12 the probabilities for these transitions are also smaller than the 10−6 K accuracy of our numerical calculations. Finally, there is the possibility that couplings such as the hyperfine interaction also yield level
repulsions with energy-level splittings that are significantly
larger than those generated by the DM interaction. We leave
the study of these interesting topics for future research.
VI. DISCUSSION

We have studied the dependence of the energy-level diagrams, with level repulsions due to the DzyaloshinskiiMoriya interaction, on the direction of the applied magnetic
field. We found that the dependence on the direction of themagnetic field seems to be generic, at least if the system has
C3 symmetry. Our numerical data suggest that the three-spin
model reproduces the main features of the 15-spin model of
V15. The presence of the Dzyaloshinskii-Moriya interaction
allows for adiabatic changes of the magnetization but, ac-

064401-8

PHYSICAL REVIEW B 70, 064401 (2004)

DZYALOSHINSKII-MORIYA INTERACTIONS AND…

cording to our calculations, the value of the resonant field forthe 01 / 2 , −1 / 2& to 0 1 / 2 , 1 / 2& transition changes with the direction of the magnetic field. The Dzyaloshinskii-Moriya interaction not only lifts the degeneracy but, depending on the
direction of the field with respect to the symmetry axis, also
shifts the resonant point away from h = 0.
The butterfly hysteresis loop observed in time-resolved
magnetization measurements has been interpreted in terms ofa combination of a Landau-Zener-Stückelberg transition atzero field and spin-phonon coupling.16,22 Our results showthat unless the field is applied in a special direction (x or y
direction in this case), the adiabatic magnetization process is
no longer symmetric with respect to the field. The dependence on the direction of the field should lead to observablechanges in the hysteresis loops. So far, only weak directional
dependence has been reported in experiments.43 Therefore,
it seems that it is necessary to explore other mechanismsthat yield energy-level repulsions such as hyperfine
interactions.52

field is parallel to the x axis, the expressions themselves arerather long and not very illuminating. Therefore they are notgiven here.
For a magnetic field parallel to the z axis, a straightforward calculation shows that
%H03/2,3/2& = −H2%03/2,3/2& =H%03/2,1/2& = −APPENDIX A: DIAGONALIZATION OF MODEL (1)Here we collect some analytical results of the solution of
the eigenvalue problem of Hamiltonian (1). We consider only
the case of a magnetic field that is parallel to the z axis 'h= !0 , 0 , hz"(. For a DM vector satisfying the conditions (2)and (3) the eight eigenvalues are given by-

E0,1 =− -3Dz + 4hz ± 9D x + 9D y + !-3Dz − 3J + 2hz",
4

E2,3 =− -3Dz − 4hz ± 9D x 2 + 9D y 2 + !-3Dz − 3J − 2hz"2
,
4

2

2

-

E4,5 =

E6,7 =-3Dz ± -3D x 2 + 3D y 2 + !-3Dz + 3J − 2hz"24-3Dz ± -3D x 2 + 3D y 2 + !-3Dz + 3J + 2hz"24z
Dz = D1,2.

!3J + 2hz"2 + 3!D2x + D2y "%03/2,1/2&
16−

3!Dx + iDy"Dz%0ā&,
8!A2"

where 0ā& denotes the state 0a& with all spins reversed and0a& =1'!1 − i-3"0 ↓ ↑ ↑ &2 3+ !1 + i-3"%0 ↑ ↓ ↑ & − 2'%0 ↑ ↑ ↓ &(.

!A3"

The expressions for H 0 3 / 2 , −3 / 2& , H2 0 3 / 2 , −3 / 2& , H 0 3 / 2 ,−1 / 2&, and H2 0 3 / 2 , −1 / 2& are obtained from Eqs. (A2) and(A3) by changing the sign of hz and Dy and replacing 0a& by
0b& =

−1

2 -3

'!1 − i-3"%0 ↑ ↓ ↓ &+ !1 + i-3"0 ↓ ↑ ↓ & − 2'%0 ↓ ↓ ↑ &(.

!A4"

Note that %a 0 b& = 0. From Eq. (A2) it follows that for theexternal field parallel to the z axis, model (1) does not allow
transitions from the state with all spins up (down) to the statewith two spins down (up). Therefore, if initially the system isin the state with all spins down, adiabatically sweeping thefield from a large negative value to a large positive value will
not yield the final state with all spins up.

,

APPENDIX B: NUMERICAL METHODS

,
!A1"

x
y
, Dy = D1,2,
Dx = D1,2

3!Dx + iDy"!-3Dz + 4hz"
%0a&,8

-3!Dx + iDy"
3J + 2hz%0ā&,%03/2,1/2& −4
4

H2%03/2,1/2& =We thank I. Chiorescu and V. Dobrovitski for illuminating
discussions. Support from the Nederlandse Stichting voor
Nationale Computer Faciliteiten (NCF) is gratefully acknowledged. The present work is partially supported by aGrant-in-Aid from the Ministry of Education, Culture,
Sports, Science and Technology, and also by the NAREGI
Nanoscience Project, Ministry of Education, Culture, Sports,
Science and Technology, Japan.

9'!J + 2hz"2 + D2x + D2y (%03/2,3/2&
16−

ACKNOWLEDGMENTS

2

3!J + 2hz"
3!Dx + iDy"
%03/2,3/2& +
%0a&,
4
4

and
Substituting thewherevalues of all model parameters, we recover the results obtained by numerical diagonalization. For hz = 0 there are four
pairs of twofold degenerate levels. Although it is possible tofind analytical expressions for the case that the magnetic

A theoretical description of quantum dynamical phenomena in the Mn12 and V15 nanomagnets requires detailed
knowledge of their energy-level schemes. Disregarding thefascinating physics of the nanomagnets, the calculation of
the eigenvalues of their model Hamiltonians is a challenging
problem in its own right. First, the (adiabatic) quantum dynamics of these systems is mainly determined by the (tiny)level repulsions. Therefore the calculation of the energy levels of these many-spin Hamiltonians has to be very accurate
in order to bridge the energy scales involved (e.g., from

064401-9

PHYSICAL REVIEW B 70, 064401 (2004)

De RAEDT et al.

500 K to# 10−19 K). Second, the level repulsions originate
from the DM interactions which mix states with different
magnetization. In principle, this prevents the use of the magnetization as a vehicle to block-diagonalize the Hamiltonianand effectively reduce the size of the matrices that have to bediagonalized. If a level repulsion involves states of significantly different magnetization (e.g., M z = −10 and M z = 10) a
perturbative calculation of the level splitting would require
going to rather high order (at least 20), a cumbersome procedure. Therefore it is of interest to explore alternative routes
to direct but accurate, brute-force diagonalization of the fullmodel Hamiltonian.
As a nontrivial set of reference data, we used the eigenvalues obtained by full diagonalization (using standard
LAPACK algorithms) of the 10 000# 10 000 matrix representing model (11).42 For one set of model parameters, such a
calculation takes about 2 h of CPU time on an Athlon
1.8 GHz, 1.5 Gbytes system. Clearly this is too slow if wewant to compute the energy-level diagram, in particular if wewant to estimate the structure of the level splittings. At theresonant fields we need the eigenvalues for many values ofh. Furthermore, in the case of V15 this calculation would takeabout 30 times longer and require about 15 Gbytes ofmemory which, for present-day computers, is too much to beof practical use.
We have tested different standard algorithms to computethe low-lying eigenvalues of large matrices. The standard
Lanczos method (including its conjugate gradient version) as
well as the power method53,54 either converge too slowly,
lack the accuracy to resolve the (nearly) degenerate eigenvalues, on sometimes even completely fail to correctly reproduce the low-lying part of the spectrum. This is not a surprise. By construction these methods work well if the groundstate is not degenerate and there is little guarantee that they
will work if there are (nearly) degenerate eigenvalues.53,54 In
particular, the Lanczos procedure suffers from numerical instabilities due to the loss of orthogonalization of the Lanczos
vectors.53,54 It seems that model Hamiltonians for the nanoscale magnets provide a class of (complex Hermitian) eigenvalue problems that are hard to solve.
Extensive tests led us to the conclusion that only the
Lanczos method with full orthogonalization53,54 and the
Chebyshev-polynomial-based projector method discussed
below can solve these rather large and difficult eigenvalue
problems with sufficient accuracy. The former is significantly
faster than the latter but using both gives extra confidence inthe results.
1. Lanczos method with full orthogonalization

In the LFO, each time a new Lanczos vector is generated
we explicitly orthogonalize (to working precision) this vector
to all, not just to the two previous, Lanczos vectors.53,54 With
some minor modifications to restart the procedure when the
Lanczos iteration terminates prematurely, after n steps thisprocedure transforms the n # n matrix H into a tridiagonal
matrix that is comparable in accuracy to the one obtained
through Householder tridiagonalization but offers noadvantages.54 In our case we are interested only in a few

low-lying eigenstates of H. Thus we can exploit the fact thatprojection onto the (numerically exact) subspace of dimension k !k * n" built by the Lanczos vectors will yield increasingly accurate estimates of the smallest (largest) eigenvalues
and corresponding eigenvectors as k increases.
In practice, to compute the M lowest energy levels, the
LFO procedure is carried out as follows.
(1) Perform a Lanczos step according to the standard procedure.
(2) Use the modified Gramm-Schmidt procedure to orthogonalize the new Lanczos vector with respect to all previous ones.53,54(3) Compute the matrix elements of the tridiagonal matrix.
(4) At regular intervals, diagonalize the tridiagonal matrix, compute the approximate eigenvectors &i , (i
= %&i 0 H 0 &i&, and %2i = %&i 0 !H − (i"2 0 &i& for i = 1 , … , M, andcheck if all %i are smaller than a specified threshold. If so,
terminate the procedure (the exact eigenvalue Ei closest to (isatisfies (i − % ' Ei ' (i + %i". If not, continue generating new
Lanczos vectors, etc.
2. Chebyshev polynomial projector method

As an alternative to the LFO, we have used a power
method53,54 based on the matrix exponential e−tH.55 Writing
the random vector +!0" in terms of the (unknown) eigenvectors /,i1 of H, we find+!t" = e−tE0',0%,00+!0"& + e−t!E1−E0",1%,10+!0"&+ e−t!E2−E0",2%,20+!0"& + ¯ (,!B1"showing limt→-+!t" / 2+!t" 2 .,0 if %,0 0 +!0"& ! 0. In thisnaive matrix-exponential version of the power method, convergence to the lowest eigenstate is exponential in t if
E 1 ) E 0.
The case of degenerate !E0 = E1 = ¯ " or very close !E0# E1 # ¯ " eigenvalues can be solved rather easily by applying the projector to a subspace instead of a single vector, in
combination with diagonalization of etH within thissubspace.55 First we fix the dimension k of the subspace bytaking k equal to or larger than the desired number of distinct
eigenvalues. The projection parameter t should be as large as
possible but nevertheless sufficiently small so that at least thefirst k terms survive one projection step. Then we generate a
set of random initial vectors +i!0" for i = 1 , … , k and set the
projection count n to zero. We compute the k lowest eigenstates by the following algorithm.55(1) Performa
projection
step+i'!n + 1"t(= e−tH+i!nt" for i = 1 , … , k.
(2) Computethek#kmatrices.
A = %+i'!n+ 1"t( 0 etH 0 +i'!n + 1"t(& = %+i'!n + 1"t( 0 +i!nt"& and B = %+i'!n+ 1"t( 0 +i'!n + 1"t(&. Note that A is Hermitian and B is positive definite.
(3) Determine the unitary transformation U that solves
the k # k generalized eigenvalue problem Ax = /Bx. Recall
that k is small.
(4) Compute +!i '!n + 1"t( = $kj=1 Ui,j+ j '!n + 1" t( for i= 1 , … , k.

064401-10

PHYSICAL REVIEW B 70, 064401 (2004)

DZYALOSHINSKII-MORIYA INTERACTIONS AND…

(5) Set +i'!n + 1"t( = +i!'!n + 1"t( for i = 1 , … , k.
(6) Compute (i = %+i'!n + 1"t( 0 H 0 +i'!n + 1"t(& and checkif %2i = %+i'!n + 1"t( 0 !H − (i"2 0 +i'!n + 1"t(& is smaller than aspecified threshold for i = 1 , … , k. If yes, terminate the calculation. If no, increase n by 1 and repeat the procedure.
We calculate e−tH+ by using the Chebyshev polynomial
expansion method.56–60 First we compute an upper bound R
of the spectral radius of H (i.e., 2H 2 ' R) by repeatedly using
the triangle inequality.60 From this point on we use the “normalized” matrix H̃ = !2H / R − 1" / 2. The eigenvalues of theHermitian matrix H̃ are real and lie in the interval'−1 , 1(.53,54 Expanding the initial value +!0" in the (unknown) eigenvectors , j of H̃ !or H" we find3+!t" = I0!z"I + 2-

4$ Im!z"Tm!H̃" +!0".

m=1

!B3"

Here, I is the identity matrix and Tm!H̃" is the matrix-valued
Chebyshev polynomial defined by the recursion relations
T0!H̃"+!0" = +!0",

T1!H̃"+!0" = H̃+!0",

!B4"

and
Tm+1!H̃"+!0" = 2H̃Tm!H̃"+!0" − Tm−1!H̃"+!0", !B5"

where z = −tR. We find the Chebyshev polynomial expansionof +!t" by computing the Fourier coefficients of the function
ezcos 0.61 Alternatively, since −1 ' Ẽ j ' 1, we can use the ex˜pansion ezE j = I0!z" + 2$m=1
Im!z"Tm!Ẽ j" where Im!z" is themodified Bessel function of integer order m (Ref. 61) towrite Eq. (B2) as

for m 1 1. In practice we will sum only contributions withm ' M where M is chosen such that, for allm ) M , 0 Im!z" / I0!z"0 is zero to machine precision. Then
it is not difficult to show that 2e−tH / I0!z" − I
M'Im!z" / I0!z"(Tm!H̃"2 is zero to machine precision too− 2$m=1[instead of e−tH we can equally well use e−tH / I0!z" as the
projector].
Using the downward recursion relation of the modified
Bessel functions, we can compute K Bessel functions to machine precision using of the order of K arithmetic
operations.61,62 A calculation of the first 20 000 modified
Bessel functions takes less than 1 s on a Pentium III
600 MHz mobile processor, using 14–15 digit arithmetic.
Hence this part of a calculation is a negligible fraction of thetotal computational work for solving the eigenvalue problem.
Performing one projection step with e−tH amounts to repeatedly using recursion (B5) to obtain T̃m!B"+!0" for k= 2 , … , M, multiply the elements of this vector by Im!z", andadd all contributions.

*Electronic

11 B. Barbara, L. Thomas, F. Lionti, A. Sulpice, and A. Caneschi, J.

˜+!t" = e−tH+!0" = ezH+!0" =$j ezE , j%, j0+!0"&,
˜

j

!B2"

address:deraedt@phys.rug.nl;
URL:
http://
www.compphys.org†Electronic address: miya@spin.t.u-tokyo.ac.jp‡
Electronic address: kristel@phys.rug.nl§
Electronic address: machida@spin.t.u-tokyo.ac.jp
1
Quantum Tunneling of Magnetization, Vol. 301 of NATO Advanced Study Institute, Series E: Applied Sciences, edited by L.
Gunther and B. Barbara (Kluwer, Dordrecht, 1995).
2 A. Caneschi, D. Gatteschi, R. Sessoli, A. Barra, L. C. Brunel, and
M. Guillot, J. Am. Chem. Soc. 113, 5873 (1991).
3 R. Sessoli, H.-L. Tsai, A. R. Shake, S. Wang, J. B. Vincent, K.
Folting, D. Gatteschi, G. Christou, and D. N. Hendrickson, J.
Am. Chem. Soc. 115, 1804 (1993).
4 D. Gatteschi, L. Pardi, A. L. Barra, and A. Müller, Mol. Eng. 3,
157 (1991).
5 G. Levine and J. Howard, Phys. Rev. Lett. 75, 4142 (1995).
6 J. R. Friedman, M. P. Sarachik, J. Tejada, and R. Ziolo, Phys.
Rev. Lett. 76, 3830 (1996).
7
L. Thomas, F. Lionti, R. Ballou, D. Gattesehi, R. Sessoli, and B.
Barbara, Nature (London) 383, 145 (1996).
8
C. Sangregorio, T. Ohm, C. Paulsen, R. Sessoli, and D. Gatteschi,
Phys. Rev. Lett. 78, 4645 (1997).
9 W. Wernsdorfer and R. Sessoli, Science 284, 133 (1999).
10 W. Wernsdorfer, T. Ohm, C. Sangregorio, R. Sessoli, D. Mailly,
and C. Paulsen, Phys. Rev. Lett. 82, 3903 (1999).

Magn. Magn. Mater. 177, 1324 (1998).
12 J. A. A. J. Perenboom, J. S. Brooks, S. Hill, T. Hathaway, and N.

S. Dalal, Phys. Rev. B 58, 330 (1998).
I. Chiorescu, W. Wernsdorfer, A. Müller, H. Bögge, and B. Barbara, Phys. Rev. Lett. 84, 3454 (2000).
14 T. Pohjola and H. Schoeller, Phys. Rev. B 62, 15 026 (2000).
15 Y. Zhong, M. P. Sarachik, J. Yoo, and D. N. Hendrickson, Phys.
Rev. B 62, R9256 (2000).
16 I. Chiorescu, W. Wernsdorfer, A. Müller, H. Bögge, and B. Barbara, J. Magn. Magn. Mater. 221, 103 (2000).
17 I. Chiorescu, R. Giraud, A. G. M. Jansen, A. Caneschi, and B.
Barbara, Phys. Rev. Lett. 85, 4807 (2000).
18 D. W. Boukhvalov, V. V. Dobrovitski, M. I. Katsnelson, A. I.
Lichtenstein, B. N. Harmon, and P. Kögerler, J. Appl. Phys. 93,
7082 (2003).
19 D. W. Boukhvalov, A. I. Lichtenstein, V. V. Dobrovitski, M. I.
Katsnelson, B. N. Harmon, V. V. Mazurenko, and V. I. Anisimov, Phys. Rev. B 65, 184435 (2002).
20 W. Wernsdorfer, N. Allaga-Alcalde, D. N. Hendrickson, and G.
Christou, Nature (London) 416, 407 (2002).
21 A. Honecker, F. Meier, D. Loss, and B. Normand, Eur. Phys. J. B
27, 487 (2002).
22 I. Chiorescu, W. Wernsdorfer, A. Müller, S. Miyashita, and B.
Barbara, Phys. Rev. B 67, 020402 (2003).
13

064401-11

PHYSICAL REVIEW B 70, 064401 (2004)

De RAEDT et al.
S. Miyashita, J. Phys. Soc. Jpn. 64, 3207 (1995).
65, 2734 (1996).
25 V. V. Dobrovitskii and A. K. Zvezdin, Europhys. Lett. 38, 377(1997).
26 L. Gunther, Europhys. Lett. 39, 1 (1997).
27 H. De Raedt, S. Miyashita, K. Saitoh, D. García-Pablos, and N.
García, Phys. Rev. B 56, 11 761 (1997).
28 S. Miyashita and N. Nagaosa, Prog. Theor. Phys. 106, 533(2001).
29 N. P. Konstantinidis and D. Coffey, Phys. Rev. B 66, 174426(2002).
30 I. Rudra, K. Saitoh, S. Ramasesha, and S. Miyashita, cond-mat/
0212354 (unpublished).
31
T. A. Kaplan, Z. Phys. B: Condens. Matter 49, 313 (1983).
32 L. Shekhtman, O. Entin-Wohlman, and A. Aharony, Phys. Rev.
Lett. 69, 836 (1992).
33
L. Shekhtman, A. Aharony, and O. Entin-Wohlman, Phys. Rev. B
47, 174 (1993).
34 A. Zheludev, S. Maslov, I. Tsukada, I. Zaliznyak, L. P. Regnault,
T. Masuda, K. Uchinokura, R. Erwin, and G. Shirane, Phys. Rev.
Lett. 81, 5410 (1998).
35 I. E. Dzyaloshinskii, Zh. Eksp. Teor. Fiz. 32, 1547 (1957) [Sov.
Phys. JETP 5, 1259 (1957)].
36
T. Moriya, Phys. Rev. 120, 91 (1960).
37 M. I. Katsnelson, V. V. Dobrovitski, and B. N. Harmon, Phys.
Rev. B 59, 6919 (1999).
38
M. Al-Saqer, V. V. Dobrovitski, B. N. Harmon, and M. I. Katsnelson, J. Appl. Phys. 87, 6268 (2000).
39 I. Rudra, S. Ramasesha, and D. Sen, Phys. Rev. B 64, 014408(2001).
40 I. Rudra, S. Ramasesha, and D. Sen, J. Phys.: Condens. Matter
13, 11717 (2001).
41
C. Raghu, I. Rudra, D. Sen, and S. Ramasesha, Phys. Rev. B 64,
064419 (2001).
42
H. De Raedt, A. H. Hams, V. V. Dobrovitsky, M. Al-Saqer, M. I.
Katsnelson, and B. N. Harmon, J. Magn. Magn. Mater. 246, 392(2002).
43 I. Chiorescu (private communication).
23

24 S. Miyashita, J. Phys. Soc. Jpn.

44

V. V. Kostyuchenko and A. K. Zvezdin, Phys. Solid State 45, 903(2003).
45 N. Regnault, Th. Jolicoeur, R. Sessoli, D. Gatteschi, and M. Verdaguer, Phys. Rev. B 66, 054409 (2002).
46 A. L. Barra, D. Gatteschi, and R. Sessoli, Phys. Rev. B 56, 8192
(1997).
47 F. Luis, F. Bartolomé, and J. Fernández, Phys. Rev. B 57, 505(1998).
48 S. Hill, J. A. A. J. Perenboom, N. S. Dalal, T. Harthaway, T.
Stalcup, and J. S. Brooks, Phys. Rev. Lett. 80, 2453 (1998).
49
E. Del Barco, A. D. Kent, E. M. Rumberger, D. N. Hendrickson,
and G. Christou, Europhys. Lett. 60, 768 (2002).
50 D. A. Garanin, J. Phys. A 24, L61 (1991).
51 D. A. Garanin and E. M. Chudnovsky, Phys. Rev. B 59, 3671
(1997).
52 S. Miyashita, H. De Raedt, and K. Michielsen, Prog. Theor. Phys.
110, 889 (2003).
53
J. H. Wilkinson, The Algebraic Eigenvalue Problem (Clarendon
Press, Oxford, 1965).
54 G. H. Golub and C. F. Van Loan, Matrix Computations (Johns
Hopkins University Press, Baltimore, 1996).
55
H. De Raedt, Comput. Phys. Rep. 7, 1 (1987).
56
H. Tal-Ezer and R. Kosloff, J. Chem. Phys. 81, 3967 (1984).
57 T. Iitaka, S. Nomura, H. Hirayama, X. Zhao, Y. Aoyagi, and T.
Sugano, Phys. Rev. E 56, 1222 (1997).
58
C. Leforestier, R. H. Bisseling, C. Cerjan, M. D. Feit, R. Friesner,
A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D.
Meyer, N. Lipkin, O. Roncero, and R. Kosloff, J. Comput. Phys.
94, 59 (1991).
59 R. N. Silver and H. Röder, Phys. Rev. E 56, 4822 (1997).
60
V. V. Dobrovitski and H. A. De Raedt, Phys. Rev. E 67, 056702(2003).
61
Handbook of Mathematical Functions, edited by M. Abramowitz
and I. Stegun (Dover, New York, 1964).
62 W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes (Cambridge University Press, New
York, 1986).

064401-12