Statistical Physics of Quantum Systems — novel orders and dynamics —J. Phys. Soc. Jpn. Vol. 74 (2005) Suppl. pp. 107–110c
2005
The Physical Society of Japan

Temperature Dependence of ESR Intensity for the Nanoscale Molecular Magnet V15
Manabu Machida1 ∗ , Toshiaki Iitaka2 † and Seiji Miyashita3 ‡1

Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan
Computational Astrophysics Laboratory, RIKEN (The Institute of Physical and Chemical Research), 2-1 Hirosawa, Wako,
Saitama 351-0198, Japan3
Department of Physics, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033,
Japan

2

(Received December 23, 2004)

The electron spin resonance (ESR) of nanoscale molecular magnet V15 is studied. Since the
Hamiltonian of V15 has a large Hilbert space and numerical calculations of the ESR signal
evaluating the Kubo formula with exact diagonalization method is diﬃcult, we implement theformula with the help of the random vector technique and the Chebyshev polynominal expansion, which we name the double Chebyshev expansion method. We calculate the temperaturedependence of the ESR intensity of V15 and compare it with the data obtained in experiment.
As another complementary approach, we also implement the Kubo formula with the subspace
iteration method taking only important low-lying states into account. We study the ESR absorption curve below 100K by means of both methods. We ﬁnd that side peaks appear due tothe Dzyaloshinsky-Moriya interaction and these peaks grows as temperature decreases.
KEYWORDS: electron spin resonance, molecular magnets

The IV V15 molecule is the complex of formulaK6 V15As6 O42 (H2 O) · 8H2 O. Since Müller and Döring
synthesized this molecule for the ﬁrst time,1–5) V15 has
been studied intensively as one of promising nanometerscale molecular magnets. In V15 , ﬁfteen 1/2 spins of
vanadium ions are arranged almost on a sphere. The
triangle cluster with three spins is sandwiched by theupper and lower hexagons. In experiment, an adiabaticchange of the magnetization has been observed in a fast
sweeping ﬁeld.6–8) In a slow sweeping ﬁeld, an interesting magnetic plateau appears. This phenomenon, called
the phonon bottleneck eﬀect, is due to an eﬀect of contact with thermal bath,6, 7) and also theoretically analyzed from a general viewpoint of the magnetic Foehn effect.9) Since the magnetization changes smoothly at zeroﬁeld from −1/2 to 1/2 when the ﬁeld is swept adiabatically at low temperatures, the gap between the groundstate and the lowest excited state at zero ﬁeld is believed to open. The most plausible origin of the gap is the
Dzyaloshinsky-Moriya (DM) interaction.8, 10–14) Around
2.8T, the ground state magnetization changes from 1/2
to 3/2. This change also occurs smoothly, and this broadness of the change is also considered to be caused by the
DM interaction.7) However, the detail of the mechanismof this broad change is not fully understood yet.15, 16)
One possible way to determine the mechanism of adiabatic change in V15 is that we compare the numericallyobtained ESR absorption curve of V15 with the experimentally obtained one. With this motivation, in this paper, we establish numerical methods to study the ESR
absorption curve obtained from the model Hamiltonianof V15 . The ESR absorption curve is calculated by the
Kubo formula in theory.17, 18) The direct implementation∗ E-mail address: machida@iis.u-tokyo.ac.jp† E-mail address: tiitaka@riken.jp‡ E-mail address: miya@spin.phys.s.u-tokyo.ac.jp

of the Kubo formula is to diagonalize the Hamiltonian
matrix.19) However, this direct method is impossible inmost cases because the dimension of the Hilbert space isquite large. For example, it is 215 = 32768 for V15 .
In this paper, we study the ESR absorption curve of
V15 by proposing two methods: one is the double Chebyshev expansion method (DCEM) and the other is the subspace iteration method (SIM). Relying on these methods,
we reveal the temperature dependence of the intensity of
V15 , and ﬁnd side peaks due to the DM interaction whichallows transitions between excited states otherwise forbidden.
The magnetic interactions in V15 are described by thefollowing Hamiltonian15, 16)



H=−
Jij S i · S j +D ij · [S i × S j ] −H S · Si.
i,j

i,j

i

(1)
We show the interactions between spins in Fig. 1. For
Jij , we have three diﬀerent values J, J1 , and J2 (|J| >|J2 | > |J1 |) with respect to the bonds on the upperand lower hexagons. Three spins on the triangle between
two hexagons interact with the hexagons by J1 and J2 .
The interactions between the three spins are negligibly
small. Here we take J = −800K, J2 = −350K, andJ1 = −225K.14) The second term on the right-hand sidein Eq. (1) describes the DM interaction. DM vectors are
considered to exist on the two hexagons in the bonds
with J. We take the reference DM vector D 1,2 to be
y
x
z
D1,2
= D1,2= D1,2= 40K.15) The other DM vectors
on the upper hexagon are obtained by rotating D1,2 byyx
2π/3 and 4π/3, i.e., D3,4
= 14.641K, D3,4
= −54.641K,yzx
D3,4 = 40K, D5,6 = −54.641K, D5,6 = 14.641K, andz
D5,6= 40K. If we assume the D3 symmetry of V15 ,1) thelower hexagon diﬀers from the upper hexagon by rotation π/6, and the z components of the DM vectors on thelower hexagon have opposite sign from those of the DM
107

108

M. Machida, et al.

J. Phys. Soc. Jpn. Vol. 74 (2005) Supplement5

Hsc . By this procedure, we obtained e−βH/2 without diagonalization. At the third step, we obtain time
evolutions of vectors |ΦBoltz ; t = e−iHt |ΦBoltz  and|ΦM x ; t = e−iHt |ΦM x , where |ΦM x  = M x |ΦBoltz  andNspin x
M x = j=1
Sj . In the BWTDM, the time evolution
is performed by the leap frog method,23) which evolves
state |φ; t as

6
6

4

5

4

3

2

1

1
8
3

9

7

2
8

10 11 12 13 14 15|φ; t + Δt = −2iHΔt |φ; t + |φ; t − Δt .

7

(3)

J

9

Note that the condition Emax Δt
1 should be satisﬁed,
where Emax is the largest eigenvalue of the Hamiltonian.
At the fourth step, we calculate the correlation function

J2

10

J1

11

15

g(t; T )
12

=

14

=
13

=
Fig. 1. The ﬁgures on the left show interactions within eachhexagon. The interaction within the triangle is negligibly small.
The ﬁgure on the top right shows interactions between the triangle and the upper and lower hexagons.

vectors on the upper hexagon. Thus, the DM vectors onx= −54.641K,
the lower hexagon are obtained as D10,11yzx
D10,11 = −14.641K, D10,11 = −40K, D12,13= 40K,yzx
D12,13= −40K, D12,13= −40K, D14,15= 14.641K,
yz= 54.641K, and D14,15= −40K. We assume
D14,15the magnetic ﬁeld is applied parallel to the c-axis of
the molecule (z-axis) and set H S = (0, 0, HS ). We take
HS = 4T throughout the paper (see also Fig. 4(b)).
Double Chebyshev expansion method (DCEM). —This method realizes O(N ) calculation with the helpof the random vector technique20, 21) and the Chebyshev polynomial expansion of exponential operators. The
DCEM is an extension of the Boltzmann-weighted timedependent method (BWTDM) developed by Iitaka andEbisuzaki.22) We brieﬂy review the BWTDM. The procedure of the BWTDM is divided by the following ﬁve
steps. By adopting the unit in Ref. 22, we have dimensionless energy, time, and temperature. At the ﬁrst step,
we prepare a random vector |Φ. For a given basis |n
of the 
Hilbert space, this random vector is given by|Φ = Nn=1 |n ξn . Here, the dimension of the Hilbert
space is N and the statistical average of the complex
random coeﬃcients {ξn } satisﬁes ξn∗  ξn  = δn n . At thesecond step, we obtain the Boltzmann-weighted random
vector |ΦBoltz  = e−βH/2 |Φ making use of the Chebyshev polynomial expansion. Here, β(= 1/T ) is the inverse
temperature. To this end, we divide the Hamiltonian byconstant Δλ as Hsc = H/Δλ so that the largest eigenvalue of Hsc does not exceed unity. Then, we expand
e−βH/2 as follows.
e−βH/2=I0 (−βΔλ/2) T0 (Hsc )

+

2

k
max

Ik (−βΔλ/2) Tk (Hsc ),

(2)

k=1

where Ik (x) is the modiﬁed Bessel function and Tk (Hsc )is the Chebyshev polynomial, which satisﬁes Tk (Hsc ) =2Hsc Tk−1 (Hsc )−Tk−2 (Hsc ), T0 (Hsc ) = 1, and T1 (Hsc ) =M x M x (t)Tr e−βH M x eiHt M x e−iHt /Tr e−βH ΦM x ; t| M x |ΦBoltz ; t , ΦBoltz |ΦBoltz  

(4)

where the trace was replaced by the inner product ofrandom vectors. Finally, the imaginary part of the dynamical susceptibility χ (ω; T ) is obtained by the Fouriertransform of g(t; T ).
 ∞


−βω
χ (ω; T ) =Re
1−eg(t; T )e−iωt dt=


1 − e−βω Re Tmax0
2 2g(t; T )e−iωt e−η t /2 dt. (5)0Here, we introduced the Gaussian ﬁlter with variance1/η 2 . This η determines the frequency resolution. The
upper limit of the integral Tmax satisﬁes Tmax ∼ 1/η inorder to avoid the Gibbs oscillation. Also η should satisfy
HS . The average energy0<η1, ηHS , and βη 2absorption per unit time I(ω; T ) is given by
ωHR2 χ (ω; T ).
(6)
2
The DCEM is almost the same as the BWTDM. Only
the third step is diﬀerent. In the DCEM, we make useof the Chebyshev polynomial expansion not only in thesecond step obtaining e−βH/2 , but also in the third stepof time evolution. State vector |φ; t is evolved with the
Chebyshev polynomials asI(ω; T ) =|φ; t + τ =e−iτ ΔλHsc |φ; t=

J0 (τ Δλ)T0 (Hsc ) |φ; t

+

2

k
max

(−i)k Jk (τ Δλ)Tk (Hsc ) |φ; t , (7)k=1

where Jk (x) is the Bessel function. Note that time stepτ is not necessarily small. In the ESR experiment for
V15 , magnetic ﬁeld HS (∼ 1K) is usually much smaller
than the strongest coupling Jmax (= |J|) between spins(∼ 103 K). That is why the frequency of precession of thespins is rather small. This means that we need to evolvestate vectors long time but do not need ﬁne resolutionof time step in order to detect small frequencies of spin
precession. Since it is possible for time step τ to havelarger value than time step Δt, the DCEM is more efﬁcient than the BWTDM for low magnetic ﬁelds. Table

M. Machida, et al.

J. Phys. Soc. Jpn. Vol. 74 (2005) Supplement
Table I. Typical computation times of two methods at various
magnetic ﬁelds.
1000(T)
100(T)
10(T)Chebyshev
126(min)
187(min)
430(min)
Leap-frog
11(min)165(min)1326(min)Intensity [T]

101
10
10

V15

0

DCEM
SIM(8)

10-2
10-3
10

levels because states of the system are conﬁned near the
ground state. The SIM is a numerical method that takes
a small subspace of the total Hilbert space which concerns only low-lying states.25, 26) Within this subspace,
the ESR absorption curve is obtained by the explicit
formulation of the Kubo formula. That is, by direct diagonalization of the small Hamiltonian. We obtain the
eigenvalues and eigenvectors of the subspace and obtain
the imaginary part of the susceptibility χ (ω);19)π   −βEmeχ (ω) =− e−βEn | ψm | M x |ψn  |2
Z m,n×δ(ω − (En − Em )),

-1

-4

1

10
100
1000
Temperature [K]10000

Fig. 2. The temperature dependence of the intensity of V15 isshown. Solid circles denote data by the DCEM and open squares
denote data by the SIM with the lowest eight levels. Dashed lines
denote, from the bottom, intensities of independent 1, 2, 3, and15 spins, respectively.

I shows typical computation times with the DCEM and
the BWTDM. Using asymptotic behavior of the Bessel
function, the ratio of two computation times is evaluated
as


b
aHS
[Chebyshev]
∼
ln
,
(8)
[Leap-frog]
Jmax
HS /Jmax
where a and b are constants.
We ﬁrst study
 the* ∞temperature dependence of the intensity I(T ) = 0 I(ω; T )dω using the DCEM. In
Fig. 2, we show the temperature dependence of I(T ).
The intensity I(T ) obtained by the DCEM is denotedby solid circles. The open squares in Fig. 2 denote dataobtained by the subspace iteration method, which wewill explain below. At very high temperatures, each spinin V15 acts as an isolated spin. Therefore, the intensity
should be given by the multiplication of I1 (T ) by 15,
where I1 (T ) = π8 HR2 HS tanh (βHS /2) is the intensity of
an isolated spin. In Fig. 2, 15×I1 (T ), 3×I1 (T ), 2×I1 (T ),and 1 × I1 (T ) are also shown by dash-dotted line, dotted line, short-dashed line, and dashed line, respectively.
As temperature decreases, the eﬀective number of spinschanges from 15 to 3. Note that the ground state magnetization is 3/2 in the present case of HS = 4T (seealso the energy diagram in Fig. 4(b)). In experiment,
the intensity starts deviating the curve for three spins ataround 200K.24) Here we conﬁrm that the intensity actually follows the curve for independent 15 spins at hightemperatures.
Subspace iteration method (SIM). — We explain another method obtaining the ESR absorption curve oflarge molecules. At low temperatures, the ESR absorption occurs only from transitions among low-lying energy

109

(9)

where Z is the partition function.
The calculation in the DCEM is numerically exact, but
it generally requires longer computation time than thatof the SIM. On the other hand, the calculation in the
SIM is carried out shorter time for a moderate value ofÑ , which is enough to study the property at low temperatures. However, it ignores higher states completely.
The two methods are complementary to each other.
In Fig. 2, the intensities obtained by the SIM with thelowest eight levels (Ñ = 8) are plotted by open squares.
Those are consistent with the solid circles obtained by
the DCEM at temperatures lower than 200K. Although
the data of the DCEM start deviating from the line ofthree spins at around 200K, the data of the SIM stayat the line, which reﬂects the fact that the SIM considers only the lowest eight levels. Here we conﬁrm that the
SIM can reproduce the data of the DCEM at low temperatures with a shorter computation time. Thus, in Fig. 2,
we can have more points at low temperatures.
We point out the following two points, which will be
discussed in more detail elsewhere.27) First, the dataslightly deviate from the dotted line at low temperatures as seen in the experiment by Sakon et al.28) For thepresent ﬁeld (HS = 4T), i.e., the ground state magnetization is 3/2, the ratio I(T )/I1 (T ) takes values smaller
than 3 at low temperatures due to the DM interaction.
Second, in the experiment by Ajiro, et al.,24) the intensity
follows the line of an isolated 1/2 spin (I(T )/I1 (T ) = 1).
This is because the ground state magnetization is 1/2 inthe ﬁeld weaker than 2.8T.
Finally, let us study the absorption curve I(ω; T ). Figure 3 compares I(ω; T = 32K) and I(ω; T = 64K) obtained by the DCEM. Note that the widths of the peaks
are solely due to the ﬁnite frequency resolution of the
DCEM given by η. We see that the peaks grow as temperature decreases. In the upper panel in Fig. 3, intensities
with and without the DM interaction are compared. We
conclude that the side peaks near the main peak around
110GHz( 4T) appear due to the DM interaction.
Figures 4(a) and 4(b) show results obtained by theSIM. In Fig. 4(a), the ESR absorptions by the lowest
eight levels (Ñ = 8) are calculated at 32K. The height ofpulses shows the value of coeﬃcients of delta functionsof the absorptions. The peak around 110GHz consists ofthree resonant peaks, whereas each side peak consists ofsingle resonant peak. Figure 4(b) shows the lowest eightenergy levels of V15 as a function of ﬁeld HS . Each la-

110

M. Machida, et al.

J. Phys. Soc. Jpn. Vol. 74 (2005) Supplement

3
32K

2
1
0
80

90

100

110
ω[GHz]

120

130

140

100

110
ω[GHz]

120

130

140

3
64K

2
1
0
80

90

Fig. 3. The upper and lower panels show the ESR absorption
curves of V15 at 4T and at 32K and 64K, respectively. In theupper panel, the dashed curve shows the absorption curve obtained from the Hamiltonian without the DM interaction.

(a)
c
10

10

d

b

5
0
109.4

5

109.6 109.8
ω [GHz]

e

a
0
70

80

-3645

90

100 110
ω [GHz]

120

140

(b)

-3650
Energy [K]

130

d

-3655
a
-3660

c e
b

-3665
0

1

2

3
Hs [T]

4

5

6

Fig. 4. (a) The ESR absorptions calculated by the SIM at 4T
and 32K are shown. The height of the pulse at each resonant frequency expresses the value of the coeﬃcient of the delta functionof the absorption. The inset shows the magniﬁed ﬁgure around
110GHz. (b) The lowest eight energy levels of V15 . The verticaldashed line at HS = 4T shows the magnetic ﬁeld applied to V15 .
The arrows denote transitions corresponding peaks labeled by athrough e in (a).

beled peak in Fig. 4(a) comes from the corresponding la-beled transition between energy levels shown in Fig. 4(b).
Thus, the combination of the DCEM and the SIM provides a powerful method to elucidate the ESR of large
molecular magnets such as V15 .
This work is supported by the Grant-in-Aid from
the Ministry of Education, Culture, Sports, Science andTechnology, and also by NAREGI Nanoscience Project,
Ministry of Education, Culture, Sports, Science andTechnology, Japan. The simulations were partially carried out by using the computational facilities of the Super Computer Center of Institute for Solid State Physics,
the University of Tokyo, and Advanced Center for Computing and Communication, RIKEN (The Institute ofPhysical and Chemical Research).
1) A. Müller and J. Döring: Angew. Chem. Int. Ed. Engl. 27 (1988)
1721.
2) A. Müller, J. Döring, and M. Penk: Z. Anorg. Allg. Chem. 595(1991) 251.
3) D. Gatteschi, L. Pardi, A. L. Barra, A. Müller and J. Döring:
Nature 354 (1991) 465.
4) A. L. Barra, D. Gatteschi, L. Pardi, A. Müller and J. Döring:
J. Am. Chem. Soc. 114 (1992) 8509.
5) D. Gatteschi, L. Pardi, A. L. Barra, and A. Müller: Mol. Eng.
3 (1993) 157.
6) I. Chiorescu, W. Wernsdorfer, A. Müller, H. Bögge, and B
Barbara: Phys. Rev. Lett. 84 (2000) 3454.
7) I. Chiorescu, W. Wernsdorfer, A. Müller, H. Bögge, and B
Barbara: J. Mag. Mag. Mat. 221 (2000) 103.
8) I. Chiorescu, W. Wernsdorfer, A. Müller, S. Miyashita, and B
Barbaraet: Phys. Rev. B 67 (2003) 020402.
9) K. Saito and S. Miyashita: J. Phys. Soc. Jpn. 70 (2001) 3385.
10) I. Dzyaloshinsky: J. Phys. Chem. Solids 4 (1958) 241.
11) T. Moriya: Phys. Rev. Lett. 4 (1960) 228.
12) T. Moriya: Phys. Rev. 120 (1960) 91.
13) S. Miyashita and N. Nagaosa: Prog. Theor. Phys. 106 (2001)
533.
14) N. P. Konstantinidis and D. Coﬀey: Phys. Rev. B 66 (2002)
174426.
15) H. De Raedt, S. Miyashita, and K. Michielsen: Phys. Stat. Sol.
B 241 (2004) 1180.
16) H. De Raedt, S. Miyashita, K. Michielsen, and M. Machida:Phys. Rev. B 70 (2004) 064401.
17) R. Kubo and K. Tomita: J. Phys. Soc. Jpn. 9 (1954) 888.
18) R. Kubo: J. Phys. Soc. Jpn. 12 (1957) 570.
19) S. Miyashita, T. Yoshino and A. Ogasahara: J. Phys. Soc. Jpn.
68 (1999) 655.
20) A. Hams and H. De Raedt: Phys. Rev. E 62 (2000) 4365.
21) T. Iitaka and T. Ebisuzaki: Phys. Rev. E 69 (2004) 057701.
22) T. Iitaka and T. Ebisuzaki: Phys. Rev. lett. 90 (2003) 047203.
23) T. Iitaka: Phys. Rev. E 49 (1994) 4684.
24) Y. Ajiro, Y. Inagaki, H. Itoh, T. Asano, Y. Narumi, K. Kindo,
T. Sakon, H. Nojiri, M. Motokawa, A. Cornia, D. Gatteschi,
A. Müller, and B. Barbara: Physica B 329-333 (2003) 1138.
25) F. Chatelin: Valeurs propres de matrices (Eigenvalues of matrices) (Masson, Paris, 1988) [in French].
26) A. Mitsutake, T. Iitaka, and Y. Okamoto: Comp. Phys. Comm.
96 (1996) 217.
27) M. Machida, T. Iitaka, and S. Miyashita: in preparation.
28) T. Sakon, K. Koyama, M. Motokawa, Y. Ajiro, A. Müller, and
B. Barbara: Physica B 346-347 (2004) 206.