Typeset with jpsj2.cls <ver.1.2> Full Paper

Fully Relativistic Full-Potential Calculations of Magnetic Moments in

Uranium Monochalcogenides with the Dirac Current

Shugo Suzuki and Hidehisa Ohta

Institute of Materials Science, University of Tsukuba, Tsukuba 305-8573

(Received April 26, 2010)

We study the orbital, spin, and total magnetic moments in uranium mono chalcogenides,

UX where X=S, Se, and Te, using the fully relativistic full-potential calculations based on

the spin density functional theory. In particular, the orbital magnetic moments are calculated

with the Dirac current. We employ two methods which adopt distinctly diﬀerent basis sets;

one is the fully relativistic full-potential linear-combination-of-atomic-orbitals (FFLCAO)

method and the other is the fully relativistic full-potential mixed-basis (FFMB) method.

Showing that the orbital magnetic moments calculated using the FFLCAO method and

those calculated using the FFMB method agree very well with each other, we demonstrate

that, in contrast to the conventional method, the method with the Dirac current enables us

to calculate the orbital magnetic moments even if the basis set includes basis functions with

no deﬁnite angular momenta, e.g., the plane waves in the FFMB method. Furthermore, it

is found that the orbital magnetic moments obtained in this work are larger by nearly 0.4

µ

B

than those obtained using the conventional method. This is crucial because the resultant

diﬀerences in the total magnetic moments are about 30 %. We compare the results of this

work with those of previous theoretical and experimental studies.

KEYWORDS: Dirac current, orbital magnetic moment, fully relativistic calculations, full-

potential calculations, uranium monochalcogenide, spin density functional theory

1. Introduction

For more than a half century, the properties of actinide compounds have been studied both

experimentally and theoretically.

1)

Among them, uranium monochalcogenides, UX where X =

S, Se, and Te, have been studied extensively as a typical material. In particular, their magnetic

properties have attracted much attention.

2–11)

The experimental studies have revealed that

UX are ferromagnetic at low temperatures. One remarkable feature of the ferromagnetism in

UX is that it is the orbital magnetic moments, M

orb

, that dominate in the total magnetic

moments, M

tot

, overcoming the spin magnetic moments, M

spin

.

12)

The calculation of M

orb

is subtle in contrast to the calculation of M

spin

; for the latter

quantity, we need only to integrate the spin density. So far, for calculating M

orb

, a conventional

method has been widely used.

13)

In this method, the magnetic-moment operator is deﬁned

in an appropriate way, and M

orb

are then calculated using the expectation values of the

magnetic-moment operator in the Bloch states. In actual calculations, the evaluation of the

1/12

J. Phys. Soc. Jpn. Full Paper

elements of the magnetic-moment operator requires little more than a rearrangement of the

overlap matrix. The magnetic moments of the unit cell calculated using this method can be

partitioned into the magnetic moments of the constituent atoms or atomic orbitals attributing

each basis function to the atom or atomic orbital to which the basis function belongs. A

disadvantage of this method is that one cannot use basis functions with no deﬁnite angular

momenta, e.g., plane waves.

Another method for calculating the magnetic moments is to use the current density as

described in the textbooks of electrodynamics.

14)

That is, if the current density is obtained,

one can calculate the magnetic moments by integrating the cross product between the posi-

tion vector and the current density. For UX as well as other actinide compounds, since the

relativistic eﬀects are signiﬁcant, the calculation of the current density should be p erformed

using the Dirac current because this includes all the relativistic eﬀects. In the fully relativis-

tic calculations based on the density functional theory adopting a single-particle equation of

the Kohn-Sham-Dirac type, the Dirac current is calculated simply using the Dirac matrices.

Furthermore, in contrast to the conventional method for calculating M

orb

, this method has

an advantage that the procedure can be applied even if the basis set includes basis functions

with no deﬁnite angular momenta. This is favorable because the physical quantities should be

calculated whatever the basis set is if its quality is good. However, to our knowledge, the cal-

culation of the magnetic moments in UX as well as other actinide compounds with the Dirac

current has not been reported so far. Thus, it seems interesting to compare M

orb

calculated

using the conventional method and those calculated using the method with the Dirac current.

When applying the method with the Dirac current, the following point should be noted.

Since the integral of the cross product between the position vector and the Dirac current does

not converge if the integral is performed over an inﬁnitely extended system, as is the case for

a crystalline solid, because of the position vector in the integrand. For this reason, an appro-

priate atomic partitioning scheme is needed. One natural choice is to use the Voronoi cells. In

actual calculations, since the Voronoi cells with sharp boundaries are not suitable for accurate

numerical calculations, the Voronoi cells with smooth boundaries are useful instead.

15–17)

It is worth pointing out that, strictly speaking, there is no guarantee of reproducing M

orb

if one employs the spin density functional theory (SDFT), the framework used widely so far,

in which only the electron density and the spin density are taken as basic variables. The theory

that takes the current density as an additional basic variable has been developed, known as

the current density functional theory (CDFT).

18–23)

Although, even within SDFT, M

orb

in

UX induced by spin-orbit coupling largely contribute to M

tot

, M

orb

are most likely enhanced

considerably when taking account of the exchange-correlation eﬀects due to the current density

as shown for 3d ferromagnetic metals.

24)

On the other hand, even if we restrict ourselves within

SDFT, it seems important to compare M

orb

calculated using the two diﬀerent methods, i.e.,

2/12

J. Phys. Soc. Jpn. Full Paper

the conventional method and the method with the Dirac current; since all the calculations of

M

orb

in UX reported so far have been performed within SDFT, the restriction within SDFT

at this stage may be useful for unambiguous comparison and also for step-by-step progress.

In this work, we study M

orb

, M

spin

, and M

tot

in UX using the fully relativistic full-

potential calculations. In particular, M

orb

are calculated with the integral of the cross product

between the position vector and the Dirac current. The calculations are performed with the

fully relativistic full-potential linear-combination-of-atomic-orbitals (FFLCAO) method and

the fully relativistic full-potential mixed-basis (FFMB) method, both based on SDFT within

the local spin density approximation (LSDA).

25)

With the two methods, whose basis sets are

distinctly diﬀerent from each other, we can examine the reliability of the results with respect

to the quality of basis sets. In §2, we describe the method of calculations. The results and

discussion are given in §3. Here, we compare the results of the FFLCAO calculations and those

of the FFMB calculations. We also compare the results of this work with those of previous

theoretical and experimental studies. Finally, we give the conclusions of this work in §4.

2. Method of Calculations

We begin with the following self-consistent equations:

[

cα · p + (β − I)mc

2

+ V

es

(r) + V

xc

(r) + βΣ · B

xc

(r)

]

ψ

ν

(r) = ε

ν

ψ

ν

(r) ,

(1)

ρ(r) = −e

∑

ν

f

ν

ψ

ν

(r)

†

ψ

ν

(r) ,

(2)

and

m(r) = −e

∑

ν

f

ν

ψ

ν

(r)

†

βΣψ

ν

(r) .

(3)

In the Dirac Hamiltonian in the left-hand side of eq. (1), c and m denote the speed of light and

the rest mass of an electron, respectively, and the rest energy of an electron, mc

2

, is subtracted.

Also, α and β are the Dirac matrices in the usual representation.

26)

In the self-consistent equa-

tions, the four-component spinor ψ

ν

(r) is the one-electron wave function of the νth level with

the energy eigenvalue ε

ν

and the occupation number f

ν

; for a crystalline solid, ν represents

the band index n and the wave vector k. In eq. (1), V

es

(r) is the electrostatic potential orig-

inated in the nuclear charges and the electron charge density, where the latter is denoted by

ρ(r) in eq. (2) with e being the electron charge. Also, V

xc

(r) = [V

up

xc

(r) + V

down

xc

(r)]/2 is the

eﬀective scalar potential that describes the spin-independent part of the exchange-correlation

potential and B

xc

(r) = [V

up

xc

(r) − V

down

xc

(r)]/2 e

z

is the eﬀective magnetic ﬁeld that describes

the spin-dependent part of the exchange-correlation potential, where V

up

xc

(r) and V

down

xc

(r)

represent the exchange-correlation potentials for up- and down-spin electrons, respectively,

and e

z

represents the unit vector along the z axis; B

xc

(r) is originated in the spin magnetiza-

tion density, m(r), which is calculated with Σ = I

2

⊗ σ where I

2

is the 2×2 unit matrix and

3/12

J. Phys. Soc. Jpn. Full Paper

σ are the usual 2×2 Pauli spin matrices. The electron charge density ρ(r) and the spin mag-

netization density m(r) are calculated with ψ

ν

(r) and f

ν

. The Dirac current is then obtained

with the following equation:

j(r) = −e

∑

ν

f

ν

ψ

ν

(r)

†

cαψ

ν

(r) .

(4)

It is crucial to note that j(r) consists of not only the orbital contribution but also the spin

contribution according to the Gordon decomposition.

27)

To divide m(r) and j(r) into atomic components, we use the atomic partitioning scheme

adopting the Voronoi cells with smo oth boundaries.

15–17)

In this scheme, the weight function

associated with the ath atom, w

a

(r), is introduced as follows:

∑

a

w

a

(r) = 1 ,

(5)

where

w

a

(r) = p(r − r

a

)/

∑

b

p(r − r

b

)

(6)

with p(r) being a function which typically is large for small arguments and small for large

arguments. In this work, we use p(r) = [exp(1/nr) − 1 − 1/nr]

n

with n = 5; instead of taking

the limit n → ∞, we take n = 5 for performing the numerical integration accurately. Using

w

a

(r), a function of space variables, f(r), is divided into atomic components, f

a

(r), as follows:

f(r) =

∑

a

f

a

(r) ,

(7)

where f

a

(r) = w

a

(r)f(r). The integral of f(r) over the whole solid, I, is also divided into

atomic components, I

a

, as follows:

I =

∑

a

I

a

,

(8)

where

I =

∫

f(r) dr (9)

and

I

a

=

∫

f

a

(r) dr . (10)

Thus, using the spin magnetization associated with the ath atom, m

a

(r) = w

a

(r)m(r), we

calculate the atomic spin magnetic moment, M

spin

a

:

M

spin

a

=

∫

m

a

(r) dr . (11)

Also, using the Dirac current associated with the ath atom, j

a

(r) = w

a

(r)j(r), we calculate

the atomic total magnetic moment, M

tot

a

:

M

tot

a

=

1

2c

∫

r × j

a

(r) dr . (12)

4/12

J. Phys. Soc. Jpn. Full Paper

Here, it may be worth mentioning again that j(r) consists of not only the orbital contribu-

tion but also the spin contribution. Accordingly, M

tot

a

consists of both the orbital and spin

contributions. Finally, we calculate the atomic orbital magnetic moment, M

orb

a

:

M

orb

a

= M

tot

a

− M

spin

a

.

(13)

An important point to be noted is that M

tot

a

calculated with eq. (12) is independent of the

choice of the origin of the position vector only if

∫

j

a

(r) dr = 0 . (14)

We have checked that this condition is always satisﬁed for the calculated results given in the

next section.

We here remark that the conventional method used in previous theoretical studies for

calculating M

orb

is diﬀerent from that used in this work although the method for calculating

M

spin

is the same. The formula used previously for calculating M

orb

is the following:

13)

M

orb

a

= −e Re

(

∑

p∈a

∑

q

∑

ν

f

ν

C

∗

pν

C

qν

∫

χ

p

(r)

†

βlχ

q

(r)dr

)

.

(15)

Here χ

p

(r) are the basis functions employed in the calculations and C

pν

are the coeﬃcients in

the expansion of ψ

ν

(r) with χ

p

(r). Also, l represent the angular momentum operator, r × p.

It is important to note that eq. (15) is applicable only if we can evaluate lχ

q

(r) deﬁnitely;

for example, it is impossible to calculate M

orb

with eq. (15) if the basis set includes plane

waves. On the contrary, we can use eqs. (11)-(13) for calculating M

orb

with any type of basis

function if the quality of the basis set is good.

UX crystallize in the NaCl structure exhibiting a strong magnetic anisotropy with an

easy axis in the [111] direction.

1)

The experimental lattice constants of US, USe, and UTe

are 5.489, 5.740, and 6.155

˚

A, respectively. These exp erimental lattice constants were used in

our calculations. We assumed that the magnetization axis is in the [111] direction, which was

taken as the z axis in our calculations. The basis functions adopted in the FFLCAO method

consist of the following four-component atomic orbitals: 1s, 2s, 2p, 3s, 3p, 3d, 4s, 4p, 4d, 4f,

5s, 5p, 5d, 5f, 6s, 6p, 6d, and 7s orbitals of the neutral U atom, 5f, 7s, and 7p orbitals of

the U

2+

atom, 1s, 2s, 2p, 3s, and 3p atomic orbitals of the neutral S atom, and 3s, 3p, and

3d orbitals of the S

2+

atom, 1s, 2s, 2p, 3s, 3p, 3d, 4s, and 4p atomic orbitals of the neutral

Se atom, and 4s, 4p, and 4d orbitals of the Se

2+

atom, 1s, 2s, 2p, 3s, 3p, 3d, 4s, 4p, 4d, 5s,

and 5p atomic orbitals of the neutral Te atom, and 5s, 5p, and 5d orbitals of the Te

2+

atom.

Also, the basis functions adopted in the FFMB method consist of the four-component atomic

orbitals of neutral U and X atoms used in the FFLCAO method and four-component plane

waves, which are positive-energy solutions of the Dirac equation for a free electron. In this

work, we chose the cut-oﬀ energy of the four-component plane waves to be 50 eV. This cut-oﬀ

energy corresponds to about 40, 50, and 60 four-component plane waves per each k point

5/12