Date

2010-04-18

Description

We present theoretical studies of the influence of spin orbit coupling on the

spin wave excitations of the Fe monolayer and bilayer on the W(110) surface.

The Dzyaloshinskii-Moriya interaction is active in such films, by virtue of the

absence of reflection symmetry in the plane of the film. When the magnetization

is in plane, this leads to a linear term in the spin wave dispersion relation

for propagation across the magnetization. The dispersion relation thus assumes

a form similar to that of an energy band of an electron trapped on a

semiconductor surfaces with Rashba coupling active. We also show SPEELS

response functions that illustrate the role of spin orbit coupling in such

measurements. In addition to the modifications of the dispersion relations for

spin waves, the presence of spin orbit coupling in the W substrate leads to a

substantial increase in the linewidth of the spin wave modes. The formalism we

have developed applies to a wide range of systems, and the particular system

explored in the numerical calculations provides us with an illustration of

phenomena which will be present in other ultrathin ferromagnet/substrate

combinations.

spin wave excitations of the Fe monolayer and bilayer on the W(110) surface.

The Dzyaloshinskii-Moriya interaction is active in such films, by virtue of the

absence of reflection symmetry in the plane of the film. When the magnetization

is in plane, this leads to a linear term in the spin wave dispersion relation

for propagation across the magnetization. The dispersion relation thus assumes

a form similar to that of an energy band of an electron trapped on a

semiconductor surfaces with Rashba coupling active. We also show SPEELS

response functions that illustrate the role of spin orbit coupling in such

measurements. In addition to the modifications of the dispersion relations for

spin waves, the presence of spin orbit coupling in the W substrate leads to a

substantial increase in the linewidth of the spin wave modes. The formalism we

have developed applies to a wide range of systems, and the particular system

explored in the numerical calculations provides us with an illustration of

phenomena which will be present in other ultrathin ferromagnet/substrate

combinations.

Type

Identifier

doi:10.1103/PhysRevB.82.014428

Phys. Rev. B 82, 014428 (2010)

Database

Link to record

Show preview

Hide preview

ar
X

iv :1

00 4.

30 66

v1 [

co nd

-m at.

me s-h

all ]

18 A

pr 20

10

Spin Orbit Coupling and Spin Waves in Ultrathin Ferromagnets:

The Spin Wave Rashba Effect

A. T. Costa,1 R. B. Muniz,1 S. Lounis,2 A. B. Klautau,3 and D. L. Mills2

1Instituto de F´ısica, Universidade Federal Fluminense, 24210-340 Nitero´i, RJ, Brasil.

2Department of Physics and Astronomy,

University of California Irvine, California, 92697, U. S. A.

3Departamento de Fisica, Universidade Federal do Para´, Bele´m, PA, Brazil.

Abstract

We present theoretical studies of the influence of spin orbit coupling on the spin wave excitations

of the Fe monolayer and bilayer on the W(110) surface. The Dzyaloshinskii-Moriya interaction is

active in such films, by virtue of the absence of reflection symmetry in the plane of the film. When

the magnetization is in plane, this leads to a linear term in the spin wave dispersion relation for

propagation across the magnetization. The dispersion relation thus assumes a form similar to

that of an energy band of an electron trapped on a semiconductor surfaces with Rashba coupling

active. We also show SPEELS response functions that illustrate the role of spin orbit coupling in

such measurements. In addition to the modifications of the dispersion relations for spin waves, the

presence of spin orbit coupling in the W substrate leads to a substantial increase in the linewidth

of the spin wave modes. The formalism we have developed applies to a wide range of systems,

and the particular system explored in the numerical calculations provides us with an illustration

of phenomena which will be present in other ultrathin ferromagnet/substrate combinations.

1

I. INTRODUCTION

The study of spin dynamics in ultrathin ferromagnets is of fundamental interest, since new

physics arises in these materials that has no counterpart in bulk magnetism. Examples are

provided by relaxation mechanisms evident in ferromagnetic resonance and Brillouin light

scattering studies,1–3 and also for the large wave vectors probed by spin polarized electron

loss spectroscopy (SPEELS).4 Of course, by now the remarkable impact of ultrathin film

structures on magnetic data storage is very well known, and other applications that exploit

spin dynamics in such materials are envisioned. Thus these issues are important from a

practical point of view as well as from that of fundamental physics.

Theoretical studies of the nature of spin waves in ultrathin films adsorbed on metal sub-

strates have been carried out for some years now, along with comparison with descriptions

provided with the Heisenberg model.5 In this paper, we extend the earlier theoretical treat-

ments to include the influence of spin orbit coupling on the spin wave spectrum of ultrathin

films. This extension is motivated by a most interesting discussion of the ground state of the

Mn monolayer on the W(110) surface. A nonrelativistic theoretical study of this system pre-

dicted that the ground state would be antiferromagnetic in character.7 This prediction was

confirmed by spin polarized scanning tunneling microscope studies of the system.8 However,

recent experimental STM data with a more sensitive instrument showed a more complex

ground state, wherein the ground state is in fact a spin density wave.9 One can construct

the new state by beginning with the antiferromagnet, and then superimposing on this a

long wavelength modulation on the direction of the moments on the lattice. The authors

of ref. 9 argued that the lack of reflection symmetry of the system in the plane of the film

activates the Dzyaloshinskii Moriya (DM) interaction, and the new state has its origin in this

interaction. They also presented relativistic and ab initio calculations that gave an excellent

account of the new data. The reflection symmetry is broken simply by the presence of the

substrate upon which the film is grown. This argument to us is most intriguing, since one

can then conclude that the DM interaction must be active in any ultrathin ferromagnet;

the substrate is surely always present. The DM interaction has its origin in the spin orbit

interaction, which of course is generally very weak in magnets that incorporate the 3d tran-

sition elements as the moment bearing entities. However, in the case of the Mn monolayer

on W(110) hybridization between the Mn 3d and the W 5d orbitals activates the very large

2

W spin orbit coupling, with the consequence that the strength of the DM interaction can

be substantial, as illustrated by the calculations presented in ref. 9. One may expect to see

substantial impact of the DM interaction in other ultrathin magnets grown on 5d substrates,

and possibly 4d substrates as well.

We have here another example of new physics present in ultrathin magnets that is not

encountered in the bulk form of the material from which the ultrathin structure is fabricated.

The purpose of this paper is to present our theoretical studies of spin orbit effects on spin

waves and also on the dynamic susceptibility of a much studied ultrathin film/substrate

combination, the Fe monolayer and bilayer on W(110). We find striking effects. For instance,

when the magnetization is in plane, as we shall see the DM interaction introduces a term

linear in wave vector in the dispersion relation of spin waves. Thus the uniform spin wave

mode at zero wave vector acquires a finite group velocity. We find this to be in the range

of 2 × 105 cm/sec for the Fe monolayer on W(110). Furthermore, left/right asymmetries

appear in the SPEELS response functions. Thus, we shall see that spin orbit coupling has

clear effects on the spin excitations of transition metal ultrathin ferromagnets grown on 5d

substrates.

We comment briefly on the philosophy of the approach used here, and in various earlier

publications.5 Numerous authors proceed as follows. One may generate a description of the

magnetic ground state of the adsorbed films by means of an electronic structure calculation

based on density functional theory. It is then possible to calculate, within the framework of

an adiabatic approximation, effective Heisenberg exchange integrals Jij between the magnetic

moments in unit cell i and unit cell j. These may be entered into a Heisenberg Hamilto-

nian, and then spin wave dispersion relations may be calculated through use of spin wave

theory. It has been known for decades10 that in the itinerant 3d magnets, effective exchange

interactions calculated in such a manner have very long range in real space. Thus, one must

include a very large number of distant neighbors in order to obtain converged results. This

is very demanding to do with high accuracy for the very numerous distant neighbors, since

the exchange interactions become very small as one moves out into distant neighbor shells.

At a more fundamental level, as noted briefly above, discussions in earlier publications

show that in systems such as we study here, the adiabatic approximation breaks down badly,

with qualitative consequences.5 First, spin wave modes of finite wave vector have very short

lifetimes, by virtue of decay into the continuum of particle hole pairs (Stoner excitations)

3

even at the absolute zero of temperature5,11 whereas in Heisenberg model descriptions their

lifetime is infinite. In multi layer films, the earlier calculations show that as a consequence of

the short lifetime, the spectrum of spin fluctuations at large wave vectors contains a single

broad feature which disperses with wave vector in a manner similar to that of a spin wave;

this is consistent with SPEELS data on an eight layer film of Co on Cu(100).4 This picture

stands in contrast to that offered by the Heisenberg model, in which a film of N layers has

N spin wave modes for each wave vector, and each mode has infinite lifetime.

The method developed earlier, and extended here to incorporate spin orbit coupling,

takes due account of the breakdown of the adiabatic approximation and also circumvents

the need to calculate effective exchange interactions in real space out to distant neighbor

shells. We work directly in wave vector space through study of the wave vector and frequency

dependent susceptibility discussed below, denoted as χ+,−( ~Q‖,Ω; l⊥, l ′ ⊥). The imaginary part

of this object, evaluated for l⊥ = l ′ ⊥ and considered as a function of frequency Ω for fixed

wave vector ~Q‖ provides us with the frequency spectrum of spin fluctuations on layer l⊥ for

the wave vector chosen. Spin waves appear as peaks in this function, very much as they do

in SPEELS data, and in a manner very similar to that used by experimentalist we extract

a dispersion relation for spin waves by following the wave vector dependence of the peak

frequency. We never need to resort to a real space summation procedure over large number

of neighbors, coupled by very tiny exchange couplings. The spin wave exchange stiffness

can be extracted either by fitting the small wave vector limit of the dispersion relation so

determined, or alternatively by utilizing an expression derived earlier5 which once again does

not require a summation in real space.

We comment on another feature of the present study. In earlier calculations,5,11,14 as in the

present paper, an empirical tight binding description forms the basis for our description of

the electronic structure. Within this approach, referred to as a multi band Hubbard model,

we can generate the wave vector and frequency dependent susceptibility for large systems.

In the earlier papers, effective tight binding parameters were extracted from bulk electronic

structure calculations. The present studies are based on tight binding parameters obtained

directly from a RS-LMTO-ASA calculation for the Fe/W(110) system. We also obtain tight

binding parameters by fitting KKR based electronic structure calculations for the ultrathin

film/substrate combinations of interest. We find that spin waves in the Fe/W(110) system

are quite sensitive to the empirical tight binding parameters which are employed, though as

4

we shall see the various descriptions provide very similar pictures of the one electron local

density of states.

We note that Udvardi and Szunyogh12 have also discussed the influence of spin orbit

coupling on the dispersion relation of spin waves in the Fe monolayer on W(110) within

the framework of the adiabatic approach discussed above, where exchange interactions and

other magnetic parameters are calculated in real space. We shall discuss a comparison with

our results and theirs below. There are differences. Most particularly, we note that in

Fig. 3, the authors of ref. 12 provide two dispersion curves for propagation perpendicular

to the magnetization, whereas in a film such as this with one spin per unit cell there can

be only one magnon branch. Additionally and very recently, Bergmann and coworkers45

investigated within an adiabatic approach finite temperature effects on the magnon spectrum

of Fe/W(110).

In section II, we comment on our means of introducing spin orbit coupling into the theory.

The results of our calculations are summarized in Section III and concluding remarks are

found in section IV.

II. CALCULATION OF THE DYNAMIC SUSCEPTIBILITY IN THE PRESENCE

OF SPIN ORBIT COUPLING

The formalism for including spin orbit coupling effects in our description of spin dynamics

is quite involved, so in this section we confine our attention to an outline of the key steps,

and an exposition of the overall structure of the theory. Our starting point is the multi band

Hubbard model of the system that was employed in our earlier study of spin dynamics in

ultrathin ferromagnets. The starting Hamiltonian is written as5

H = ∑ ij

∑ µνσ

T µν ij c

† iµσcjνσ +

1

2

∑ µνµ′ν′

∑ iσσ′

Ui;µν,µ′ν′c † iµσc

† iνσ′ciν′σ′ciµ′σ (1)

where i and j are site indices, σ, σ′ refer to spin, and µ, ν to the tight binding orbitals, nine in

number for each site, which are included in our treatment. The Coulomb interactions operate

only within the 3d orbitals on a given lattice site. The film, within which ferromagnetism

is driven by the Coulomb interactions, sits on a semi-infinite substrate within which the

Coulomb interaction is ignored.

5

In our empirical tight binding picture, the spin orbit interaction adds a term we write as

HSO = ∑ i

∑ µν

λi

2

[ Lzµν(c

† iµ↑ciν↑ − c

† iµciν↓) + L

+ µνc

† iµ↓ciν↑ + L

− µνc

† iµ↑ciν↓

] (2)

where ~L is the angular momentum operator, λi is the local spin-orbit coupling constant,

L± = Lx± iLy and Lαµν = 〈µ|L α|ν〉. We assume that the spin orbit interaction, present both

within the ferromagnetic film and the substrate, operates only within the 3d atomic orbitals.

A convenient tabulation of matrix elements of the orbital angular momentum operators is

found in ref. 13.

Information on the spin waves follows from the study of the spectral density of the trans-

verse dynamic susceptibility χ+,−( ~Q‖,Ω; l⊥, l ′ ⊥) as discussed above. From the text around

Eq. (1) of ref. 14, we see that this function describes the amplitude of the transverse spin

motion (the expectation value of the spin operator S+ in the layer labeled l⊥) to a fictitious

transverse magnetic field of frequency Ω and wave vector ~Q‖ parallel to the film surface that

is applied to layer l′⊥ of the sample. The spectral density, given by Im{χ+,−( ~Q‖,Ω; l⊥, l

′ ⊥)},

when multiplied by the Bose Einstein function n(Ω) = [exp(βΩ)−1]−1 is also the amplitude

of thermal spin fluctuations of wave vector ~Q‖ and frequency Ω in layer l⊥. We obtain in-

formation regarding the character (frequency, linewidth, and amplitude in layer l⊥) of spin

waves from the study of this function, as discussed earlier.5

Our previous analyses are based on the study of the dynamic susceptibility just described

through use of the random phase approximation (RPA) of many body theory. The Feynman

diagrams included in this method are the same as those incorporated into time dependent

density functional theory, though use of our Hubbard model allows us to solve the result-

ing equation easily once the very large array of irreducible particle hole propagators are

generated numerically.

Our task in the present paper is to extend the RPA treatment to incorporate spin orbit

coupling. The extension is non trivial. The quantity of interest, referred to in abbreviated

notation as χ+,−, may be expressed as a commutator of the spin operators S + and S− whose

precise definition is given earlier.5,12 With spin orbit coupling ignored, the RPA decoupling

procedure leads to a closed equation for χ+,−. When the RPA decoupling is carried out in

its presence, we are led to a sequence of four coupled equations which include new objects

we may refer to as χ−,−, χ↑,− and χ↓,−. The number of irreducible particle hole propagators

that must be computed likewise is increased by a factor of four. For a very simple version

6

of a one band Hubbard model, and for a very different purpose, Fulde and Luther carried

out an equivalent procedure many years ago15. In what follows, we provide a summary of

key steps along with expressions for the final set of equations.

To generate the equation of motion, we need the commutator of the operator S+µν(l, l ′) =

c † lµ↑cl′ν↓ with the Hamiltonian. One finds

[S+µν(l, l ′), HSO] =

1

2

∑ η

{λl′L + νηc

† lµ↑cl′η↑ − λlL

+ ηµc

† lη↓cl′ν↓ + λl′L

z νηc

† lµ↑cl′η↓ − λlL

z ηµc

† lη↑cl′ν↓}. (3)

The last two terms on the right hand side of Eq. 3 lead to terms in the equation of motion

which involve χ+,− whereas the first two terms couple us to the entities χ↑,− and χ↓,−.

When we write down the commutator of these new correlation functions with the spin orbit

Hamiltonian, we are led to terms which couple into the function χ−,− which is formed from

the commutator of two S− operators. In the absence of spin orbit coupling, a consequence of

spin rotation invariance of the Hamiltonian is that the three new functions just encountered

vanish. But they do not in its presence, and they must be incorporated into the analysis.

One then introduces the influence of the Coulomb interaction into the equation of motion,

and carries out an RPA decoupling of the resulting terms. The analysis is very lengthy, so

here we summarize only the structure that results from this procedure. Definitions of the

various quantities that enter are given in the Appendix. We express the equations of motion

in terms of a 4 × 4 matrix structure, where in schematic notation we let χ(1) = χ+,−,

χ(2) = χ↑,−, χ (3) = χ↓,− and χ

(4) = χ−,−. The four coupled equations then have the form

Ωχ(s) = A(s) + ∑ s′

(Bss ′

+ B˜ss ′

)χ(s ′) (4)

Each quantity in Eq. 4 has attached to it four orbital indices, and four site indices. To be

explicit, χ(2) = χ↑,− which enters Eq. (4) is formed from the commutator of the operator

c † lµ↑cl′ν↑ with c

† mµ′↓cm′ν′↑ and in full we denote this quantity as χ

(2) µν;µ′ν′(ll

′;mm′). The site

indices label the planes in the film, and we suppress reference to Ω and ~Q‖. The products

on the right hand side of Eq. 4 are matrix multiplications that involve these various indices.

For instance, the object ∑

s′ B ss′χ(s

′) is labeled by four orbital and four site indices so

[Bss ′

χ(s ′)]µν,µ′ν′(ll

′;mm′) = ∑ γδ

∑ nn′

Bss ′

µν,γδ(ll ′;nn′)χ

(s′) γδ,µ′ν′(nn

′;mm′). (5)

One proceeds by writing Eq. 4 in terms of the dynamic susceptibilities that characterize

the non-interacting system. These, referred to also as the irreducible particle hole propaga-

tors, are generated by evaluating the commutators which enter into the definition of χ(s) in

7

the non interacting ground state. These objects, denoted by χ(0s) obey a structure similar

to Eq. 4,

Ωχ(0s) = A(s) + ∑ s′

Bss ′

χ(s ′). (6)

It is then possible to relate χ(s) to χ(0s) through the relation, using four vector notation,

~χ(Ω) = ~χ(0)(Ω) + (Ω− B)−1B˜~χ(Ω). (7)

The matrix structure Γ ≡ (Ω−B)−1 may be generated from the definition of B, which may

be obtained from the equation of motion of the non-interacting susceptibility, Eq. 6. Then

B˜ follows from the equation of motion of the full susceptibility, as generated in the RPA.

One may solve Eq. 7

~χ(Ω) = [I − (Ω−B)−1B¯]−1~χ(0)(Ω), (8)

so our basic task is to compute the non interacting susceptibility matrix ~χ(0) and then carry

out the matrix inversion operation displayed in Eq. 8. For this we require the single particle

Greens functions (SPGFs) associated with our approach.

To generate the SPGFs, we set up an effective single particle Hamiltonian Hsp by intro-

ducing a mean field approximation for the Coulomb interaction. The general structure of

the single particle Hamiltonian is

Hsp = ∑ ij

∑ µνσ

T˜ µνσ ij c

† iµσcjνσ +

∑ i

∑ µν

{α∗i;µνc † iµ↓ciν↑ + αi;µνc

† iµ↑ciν↓} (9)

where the effective hopping integral T˜ µνσij contains the spin diagonal portion of the spin orbit

interaction, along with the mean field contributions from the Coulomb interaction. The form

we use for the latter is stated below. The coefficients in the spin flip terms are given by

αi;µν = λiL − µν −

∑ ηγ

Ui;ηµ,νγ〈c † iη↓ciγ↑〉. (10)

We then have the eigenvalue equation that generates the single particle eigenvalues and

eigenfunctions in the form Hsp|φs〉 = Es|φs〉; we can write this in the explicit form

∑ l

∑ ησ′

[ δσσ′ T˜

µησ′

il + δil(δσ′↓δσ↑α ∗ l;µη + δσ′↑δσ↓αl;µη)

] 〈lησ′|φs〉 = Es〈iµσ|φs〉. (11)

The single particle Greens function may be expressed in terms of the quantities that enter

Eq. 11. We have for this object the definition

Giµσ;jνσ′(t) = −iθ(t)〈{ciµσ(t), c † jνσ′(0)}〉 (12)

8

and one has the representation

Giµσ;jνσ′(Ω) = ∑ s

〈iµσ|φs〉〈φs|jνσ ′〉

Ω−Es + iη . (13)

These functions may be constructed directly from their equations of motion, which read,

after Fourier transforming with respect to time,

− ∑ l

∑ ησ′′

[ δσσ′′ T˜

µησ′′

il + δil(δσ′′↓δσ↑α ∗ l;µη + δσ′′↑δσ↓αl;µη)

] Glησ′′;jνσ′ + ΩGiµσ;jνσ′ = δσσ′δµνδij.

(14)

For the case where the substrate is semi infinite, our means of generating a numerical solution

to the hierarchy of equations stated in Eq. 14 has been discussed earlier. What remains is

to describe how the Coulomb interaction enters the effective hopping integrals T˜ µνσij that

appear in Eq. 9, Eq. 11 and Eq. 14.

There are, of course, a large number of Coulomb matrix elements in the original Hamil-

tonian, even if the Coulomb interactions are confined to within the 3d shell. Through the

use of group theory,17 the complete set of Coulomb matrix elements may be expressed in

terms of three parameters. These are given in Table I of the first cited paper in ref. 5.

In subsequent work, we have found that a much simpler structure18 nicely reproduces re-

sults obtained with the full three parameter form. We use the simpler one parameter form

here, for which Ui;µν,µ′ν′ = Uiδµν′δµ′ν . Then in the mean field approximation, the Coulomb

contribution to the single particle Hamiltonian assumes the form

H(C)sp = − ∑ i

Uimi

2

∑ µ

(c†iµ↑ciµ↑ − c † iµ↓ciµ↓) (15)

Here mi is the magnitude of the moment on site i. The Coulomb interactions Ui are non zero

only within the ultrathin ferromagnet, and the moments mi, determined self consistently,

vary from layer to layer when we consider multi layer ferromagnetic films.

It should be noted that when the Ansatz just described is employed in Eq. 10, the term

from the Coulomb interaction on the right hand side becomes proportional to the transverse

component of the moment located on site i and this vanishes identically. Thus, despite the

complexity introduced by the spin orbit coupling, when the simple one parameter Ansatz for

the Coulomb matrix elements is employed, one needs no parameters beyond the moment on

each layer in the self consistent loops that describe the ground state. In the present context,

this is an extraordinarily large savings in computational labor, and this will allow us to

9

address very large systems in the future. It is the case that certain off diagonal elements

such as 〈c†mµ′↓cl′ν↑〉 appear in the quantities defined in the Appendix. Notice, for example,

the expressions in Eqs. A.1. Once the ground state single particle Greens functions are

determined, such expectation values are readily computed.

III. RESULTS AND DISCUSSION

In earlier studies of Fe layers on W(110),19,20 as noted above, the electronic structure

was generated through use of tight binding parameters obtained from bulk electronic struc-

ture calculations. These calculations generate effective exchange interactions comparable

in magnitude to those found in the bulk transition metals,20 with the consequence that for

both monolayer Fe and bilayer Fe on W(110) the large wave vector spin waves generated

by theory are very much stiffer than found experimentally21,22 though it should be noted

that for the bilayer, the calculated value of the spin wave exchange stiffness is in excellent

accord with the data.23 Subsequent calculations which construct the spin wave dispersion

relation from adiabatic theory based on calculations of effective exchange integrals also gen-

erate spin waves for the monolayer substantially stiffer than found experimentally,12 though

they are softer than in our earlier work by a factor of two or so. We remark that it has

been suggested that the remarkably soft spin waves found experimentally may have origin

in carbon contamination of the monolayer and bilayer.20 We remark here that this can be

introduced during the SPEELS measurement. We note that the magnetic properties of Fe

monolayers grown on carbon free W(110)24 differ dramatically from those grown on surfaces

now known to be contaminated by carbon.25 In the former case, the domain walls have a

thickness of 2.15 nm,24 whereas in the latter circumstance very narrow walls with thickness

bounded from above by 0.6 nm are found.25 This suggests that the strength of the effective

exchange is very different in the two cases, with stiffer exchange in the carbon free samples.

The considerations of the previous paragraph have motivated us to carry out a series of

studies of the effective exchange in the Fe monolayer on W(110) within the framework of

three different electronic structure calculations. We find that although all three give local

density of states that are very similar, along with very similar energy bands when these are

examined, the intersite exchange interactions vary substantially. First, we have employed the

parameter set used earlier that is based on bulk electronic structures19,20 in new calculations

10

we call case A. In case B, we have employed an approach very similar to that used in ref. 12,

though in what follows our calculation of effective exchange integrals is non relativistic.

This is the Korringa Kohn Rostoker Greens Function (KKR-GF) method,26 which employs

the atomic sphere approximation and makes use of the Dyson equation G = g + gV G as

given in matrix notation. This allows us to calculate the Greens function G of an arbitrary

complex system given the perturbing potential V and the Greens function g of a reference

unperturbed system. Within the Local Spin Density Approximation (LSDA),27 We consider

a slab of five monolayers of W with the experimental lattice constant on top of which an Fe

monolayer is deposited and relaxed by -12.9%12 with respect to the W interlayer distance.

Angular momenta up to lmax = 3 were included in the Greens functions with a k mesh of

6400 points in the full two dimensional Brillouin zone. The effective exchange interactions

were calculated within the approximation of infinitesmal rotations28 that allows one to use

the magnetic force theorem. This states that the energy change due to infinitesmal rotations

in the moment directions can be calculated through the Kohn Sham eigenvalues.

Method C is the Real Space Linear-Muffin-Tin-Orbital approach as implemented, also,

in the atomic sphere approximation (RS-LMTO-ASA).29–33 Due to its linear scaling, this

method allows one to address the electronic structure of systems with a large number of

atoms for which the basic eigenvalue problem is solved in real space using the Haydock

recursion method. The Fe overlayer on the W(110) substrate was simulated by a large

bcc slab which contained ∼6800 atoms, arranged in 12 atomic planes parallel to the (110)

surface, with the experimental lattice parameter of bulk W. One empty sphere overlayer

is included, and self consistent potential parameters were obtained for the empty sphere

overlayer, the Fe monolayer, and the three W layers underneath using LSDA.34 For deeper

W layers we use bulk potential. Nine orbitals per site (the five 3d and 4 sp complex) were

used to describe the Fe valence band and the empty sphere overlayer, and for W the fully

occupied 4f orbitals were also included in the core. To evaluate the orbital moments we use

a scalar relativistic (SR) approach and include a spin orbit coupling term λ~L · ~S at each

variational step.35 In the recursion method the continued fraction has been terminated after

30 recursion levels with the Beer Pettifor terminator.36 The TB parameters so obtained are

inserted into our semi-empirical scheme and this allows us to generate the non interacting

susceptibilities which enter our full RPA description of the response of the structure.

In order to compare the electronic structures generated by the approaches just described,

11

we turn our attention to the local density of states for the majority and minority spins in

the adsorbed Fe monolayer. These are summarized in Fig. 1.

The local densities of states (LDOS) generated by the three sets of TB parameters have

approximately the same overall features, as we see from Fig. 1. The main differences appear

in the majority spin band, which overlaps the 5d states in the W substrate over a larger

energy range than the minority band. This is also true if we compare the LDOS generated

by the tight binding parameters extracted from the KKR electronic structure to the LDOS

obtained directly from the KKR calculations (red dashed line in Fig. 1b). The Fe-W hopping

parameters are indeed the least accurate portion of our parametrization scheme. In case A

we just used the Fe-Fe bulk parameters to describe the Fe-W hopping. In case B we extracted

TB parameters for Fe by fitting a KKR calculation of an unsupported Fe monolayer with a

lattice parameter matching that of the W substrate. For the Fe-W hopping we used the Fe

parameters obtained from the fitting, scaled to mimic the Fe-W distance relaxation. The

relaxation parameter was chosen to give the correct spin magnetic moment for the adsorbed

Fe monolayer. In case C all parameters were directly provided by the RS-LMTO-ASA code,

but in the DFT calculations the Fe-W distance was assumed to be equal to the distance

between W layers. Thus, the main difference between cases B and C is the treatment of the

mixing between Fe and W states and this is expected to affect more strongly the states that

occupy the same energy range.

As noted above, while the local density of states provided by the three approaches to

the electronic structure are quite similar as we see from Fig. 1 (and the same is true of the

electronic energy bands themselves if these are examined), the exchange interactions differ

substantially for the three cases. For the first, second and third neighbors we have (in meV)

42.5, 3.72 and 0.46 for model A, 28.7, -7.87 and 0.31, for model B and 11.23, -7.31 and 0.22

for model C. The authors of ref. 12 find 10.84, -3.34 and 3.64 for these exchange integrals.

We now turn to our studies of spin excitations in the Fe monolayer and Fe bilayer on

W(110) within the framework of the electronic structure generated through use of the ap-

proach in case C. We will discuss the influence of spin orbit coupling on both the trans-

verse wave vector dependent susceptibility though study of the spectral density function

A( ~Q‖,Ω; l⊥) = −Im{χ+,−( ~Q‖,Ω; l⊥)} discussed in section I. This function, for fixed wave

vector ~Q‖, when considered as a function of frequency Ω, describes the frequency spectrum

of the fluctuations of wave vector ~Q‖ of the transverse magnetic moment in layer l⊥ as noted

12

-6 -4 -2 0 2 E-EF (eV)

-2

0

2

4

LD O

S (st

ate s/e

V)

-2

0

2

4 LD

O S

(st ate

s/e V)

-2

0

2

4

LD O

S (st

ate s/e

V)

a)

b)

c)

FIG. 1: (color online) For the Fe monolayer adsorbed on W(110), we show the local density of

states in the Fe monolayer. The majority spin density of states is shown positive and the minority

spin density of states is negative. The zero of energy is at the Fermi energy. In (a), bulk electronic

structure parameters are used as in the second of the two papers cited in 5 (caseA). In (b), we

have the density of states generated by method B. The black curve is found by fitting the KKR

electronic energy band structure to tight binding parameters as described in the text, and the red

curve is calculated directly from the KKR calculation. In (c) we have the local density of states

generated by method C.

13

above. In the frequency regime where spin waves are encountered, this function is closely

related to (but not identical to) the response function probed in a SPEELS measurement.

In Fig. 2, for the Fe monolayer on W(110), we show the spectral density function cal-

culated for three values of | ~Q‖|, for propagation across the magnetization. Thus, the wave

vector is directed along the short axis in the surface. This is the direction probed in SPEELS

studies of the Fe monolayer on this surface.22 In each figure, we show three curves. The green

dashed curve is calculated with spin orbit coupling set to zero. We show only a single curve

for this case, because the spectral density is identical for the two directions of propagation

across the magnetization, + ~Q‖ and −~Q‖. When spin orbit coupling is switched on, for the

two directions just mentioned the response function is very different, as we see from the red

and black curve in the various panels. These spin wave frequencies, deduced from the peak in

the response functions as discussed in section I, differ for the two directions of propagation,

and also note that the peak intensities and linewidths differ as well. It is the absence of both

time reversal symmetry and reflection symmetry which renders + ~Q‖ and −~Q‖ inequivalent

for this direction of propagation. The system senses this breakdown of symmetry through

the spin orbit interaction. If one considers propagation parallel to the magnetization, the

asymmetries displayed in Fig. 2 are absent. The reason is that for this direction of propaga-

tion, reflection in the plane that is perpendicular to both the magnetization and the surface

is a good symmetry operation of the system, but takes + ~Q‖ into −~Q‖ thus rendering the two

directions equivalent. Recall, of course, that the magnetization is a pseudo vector in regard

to reflections. Notice how very broad the curves are for large wave vectors; the lifetime of

the spin waves is very short indeed.

As discussed in section I, we may construct a spin wave dispersion curve by plotting the

maxima in spectral density plots such as those illustrated in Fig. 2 as a function of wave

vector. We show dispersion relations constructed in this manner in Fig. 3, with spin orbit

coupling both present and absent. In Fig. 3a, and for propagation perpendicular to the

magnetization we show the dispersion curve so obtained for wave vectors throughout the

surface Brillouin zone, and in Fig. 3b we show its behavior for small wave vectors.

Let is first consider Fig, 3(a). Here the dispersion curve extends throughout the two

dimensional Brillouin zone. At the zone boundary, quite clearly the slope of the dispersion

curve does not vanish. In this direction of propagation, the nature of the point at the zone

boundary does not require the slope to vanish. What is most striking, clearly, is the anomaly

14

0 10 20 30 40 50 Ω (meV)

0

5×103

1×104

- Im

χ + ,−

(a)

0 50 100 150 200 Ω (meV)

0

500

1000

1500

- Im

χ + ,-

(b)

0 100 200 300 Ω (meV)

0

200

400

600

800

- Im

χ + ,-

(c)

FIG. 2: (color online) The spectral density functions A( ~Q‖,Ω; l⊥) evaluated in the Fe monolayer for

three values of the wave vector in the direction perpendicular to the magnetization. We have (a)

| ~Q‖| = 0.4A˚ −1 , (b) | ~Q‖| = 1.0A˚

−1 and (c) | ~Q‖| = 1.4A˚

−1 . The green curve (dashed) is calculated

with spin orbit coupling set to zero; the spectral density here is independent of the sign of ~Q‖. The

red and black curves are calculated with spin orbit coupling turned on. Now we see asymmetries

for propagation across the magnetization, with the red curve ~Q‖ directed from left to right and the

black curve from right to left.

in the vicinity of 1 A˚ −1 . This feature is evident in the calculation with spin orbit coupling

absent, and for positive values of the wave vector the feature becomes much more dramatic

when spin orbit coupling is switched on. Anomalies rather similar to those in the black

curve in Fig. 3(a) appear in the green dispersion curve found in Fig. 3 of ref. 12, though

these authors did not continue their calculation much beyond the 1 A˚ −1

regime. Our spin

waves are very much softer than theirs in this spectral region, notice. In Fig. 3 of ref. 12,

one finds two dispersion curves, one a mirror image of the second. Thus, these authors

15

-1.5 -1 -0.5 0 0.5 1 1.5

Q (A-1) 0

50

100

150

200

Ω

(m eV

) (a)

-0.4 -0.2 0 0.2 0.4

Q (A-1) 0

10

20

30

40

Ω

(m eV

)

(b)

FIG. 3: (color online) Spin wave dispersion relations constructed from peaks in the spectral den-

sity, for the Fe monolayer on W(110). The wave vector is in the direction perpendicular to the

magnetization. The red curve is constructed in the absence of spin orbit coupling, it is included in

the black curve.

display two spin wave frequencies for each wave vector. This surely is not correct. For a

structure with one atom per unit cell, there is one and only one spin wave mode for each

wave vector, though as discussed above for the structure explored here symmetry allow the

left/right asymmetry in the dispersion curve illustrated in our Fig. 3.

In Fig. 3b, again with spin orbit coupling switched on and off, we show an expanded

view of the dispersion curve for small wave vectors. With spin orbit interaction switched off,

at zero wave vector we see a zero frequency spin wave mode, as required by the Goldstone

theorem when the underlying Hamiltonian is form invariant under spin rotation. The curve

is also symmetrical, and is accurately fitted by the form Ω( ~Q‖) = 149Q 2 ‖ (meV), with the

wave vector in A˚ −1 , whereas with spin orbit coupling turned on the dispersion relation is

fitted by Ω( ~Q‖) = 3.4−11.8Q‖+143Q 2 ‖ (meV). Spin orbit coupling introduces an anisotropy

gap at Q‖ = 0, and most striking is the term linear in wave vector. This has its origin in the

Dzyaloshinskii Moriya interaction whose presence, as argued by the authors of ref. 9, has its

origin in the absence of both time reversal and inversion symmetry, for the adsorbed layer.

At long wavelengths, one may describe spin waves by classical long wavelength phe-

nomenology. The linear term in the dispersion curve has its origin in a term in the energy

density of the spin system of the form

VDM = −Γ ∫ dxdzSy(x, z)

∂Sx(x, z)

∂x (16)

16

Here Sα(x, z) is a spin density, the xz plane is parallel to the surface, and the magnetization

is parallel to the z direction.

One interesting feature of the spin wave mode whose dispersion relation is illustrated in

Fig. 3b is that at Q‖ = 0, the mode has a finite group velocity. The fit to the dispersion

curve gives this group velocity to be ∂Ω(~Q‖)

∂Q‖ ≈ 2×105cm/s , which is in the range of acoustic

phonon group velocities.

We turn now to our calculations of spin waves and the response functions for the Fe bilayer

on W(110). Let us first note that experimentally the orientation of the magnetization in

the bilayer appears to be dependent on the surface upon which the bilayer is grown. For

instance, when the bilayer is on the stepped W(110) surface, it is magnetized perpendicular

to the surface,25 a result in agreement with ab initio calculations of the anisotropy realized in

the epitaxial bilayer.38 However, in the SPEELS studies of spin excitations in the bilayer22,39

the magnetization is in plane. In our calculations, we find for model B the magnetization is

perpendicular to the surface, whereas in model C it lies in plane, along the long axis very

much as in the SPEELS experiments. The anisotropy in the bilayer is not particularly large,

on the order of 0.5 meV/Fe atom, and one sees from these results that it is a property quite

sensitive to the details of the electronic structure. The fact that model B and model C give

the two different stable orientation of the magnetization allows us to explore spin excitations

for the two different orientations of the magnetization.

We first turn our attention to the case where the magnetization lies in plane. The bilayer

has two spin wave modes, an acoustic mode for which the magnetization in the two planes

precesses in phase, and an optical mode for which they precess 180 degrees out of phase.

In Fig. 4, we show calculations of the dynamic susceptibility in the frequency range of the

acoustic mode for two values of the wave vector, Q‖ = +0.5A˚ − 1 and Q‖ = −0.5A˚

− 1. A

spin orbit induced left right asymmetry is clearly evident both in the peak frequency and

the height of the feature. Very recently, beautiful measurements of spin orbit asymmetries

in the Fe bilayer have appeared,39 and the results of our Fig. 4 are to be compared with

Fig. 3 of ref. 39. Theory and experiment are very similar, both in regard to the intensity

asymmetry and also the spin orbit induced frequency shift, though our calculated spin wave

frequencies are a little stiffer than those found experimentally.

As remarked above, in Fig. 4 we show only the acoustical spin wave mode frequency

regime. In Fig. 5, for the spectral densities in the innermost layer (upper panel) and the

17

0 50 100 Ω (meV)

0

500

1000

1500

2000

2500

3000

- Im

χ + ,−

FIG. 4: (color online) The spectral density in the innermost layer, in the acoustic spin wave regime,

for wave vectors of Q‖ = +0.5A˚ −1

(black curve) and Q‖ = −0.5A˚ −1

(red curve). Model C has

been used for the calculation. In the ground state, the magnetization lies in plane along the long

axis.

outermost layer (lower panel) we show the spectral densities for the entire spin wave regime,

including the region where the optical spin wave is found. It is clear that the spin orbit

induced frequency shifts are largest for the optical mode which, unfortunately is not observed

in the experiments.39

In Fig. 6 for a sequence of wave vectors, all chosen positive, we show a sequence of spectra

calculated for the entire frequency range so both the acoustic and optical spin wave feature

are displayed. The black curves show the spectral density of the innermost Fe layer, and

the red curves are for the outer layer. The optical spin wave mode, not evident in the data,

shows clearly in these figures. Notice that for wave vectors greater than 1A˚ −1

the acoustical

mode is localized in the outer layer and the optical mode is localized on the inner layer.

The optical mode is very much broader than the acoustical mode at large wave vectors, by

virtue of the strong coupling to the electron hole pairs in the W 5d bands.

An interesting issue is the absence of the optical mode from the SPEELS spectra reported

18

-1500

-1000

-500

0

- Im

χ + ,−

0 50 100 150 200 250 Ω (meV)

-2000

-1000

0

- Im

χ + ,−

FIG. 5: (color online) For the wave vector Q‖ = 0.5A˚ −1

we show the spectral densities in the

innermost Fe layer (upper panel) and in the outermost layer (lower panel) for the Fe bilayer on

W(110). The figure includes the optical spin wave feature. As in Fig. 4, the black curve is calculated

for Q‖ positive, and the red curves are for Q‖ negative. The calculations employ model C.

in refs. 22 and 39. We note that these spectra are taken with only two beam energies, 4 eV

and 6.75 eV. At such very low energies, the beam electron will sample both Fe layers,

so the SPEELS signal will be a coherent superposition of electron waves backscattered

from each layer; the excitation process involves coherent excitation of both layers by the

incident electron. As a consequence of the 180 degree phase difference in spin motions

associated with the two modes it is quite possible, indeed even probable, that for energies

where the acoustical mode is strong the intensity of the optical mode is weak, by virtue

of quantum interference effects in the excitation scattering amplitude. In earlier studies

of surface phonons, it is well documented that on surfaces where two surface phonons of

different polarization exist for the same wave vector, one can be silent and one active in

electron loss spectroscopy.40 It would require a full multiple scattering analysis of the spin

wave excitation process to explore this theoretically. While earlier41 calculations that address

SPEELS excitation of spin waves described by the Heisenberg model could be adapted for

19

0 2000 4000 6000

Im χ +

,−

0

1000

2000 Im

χ + ,−

0

500

1000

Im χ +

,−

0

500

1000

1500

Im χ +

,−

0 500

1000 1500

Im χ +

,−

0 100 200 300 400 500 Ω (meV)

0

500

1000

Im χ +

,−

Qy=0.4 A -1

Qy=0.6 A -1

Qy=0.8 A -1

Qy=1.0 A -1

Qy=1.2 A -1

Qy=1.4 A -1

FIG. 6: (color online) For the Fe bilayer and for several values of the wave vector (all positive), we

show the spectral density functions for the innermost layer adjacent to the substrate (black curve)

and those for the outer layer of the film. The calculations employ model C.

20

-1.5 -1 -0.5 0 0.5 1 1.5 Q (A-1)

0

100

200

300

Ω

(m eV

)

FIG. 7: (color online) For the Fe bilayer with magnetization in plane, we show the spin wave

dispersion curves calculated with spin orbit coupling (black points) and without spin orbit coupling

(red curves). Model C has been employed for these calculations.

this purpose, in principle, a problem is that at such low beam energies it is necessary to take

due account of image potential effects to obtain meaningful results.42 This is very difficult

to do without considerable information on the electron reflectivity of the surface.42 It would

be of great interest to see experimental SPEELS studies of the Fe bilayer with a wider range

of beam energies to search for the optical mode, if this were possible.

In Fig. 7, we show dispersion curves for the optical and acoustic spin wave branches for

the bilayer. The magnetization lies in plane, and one can see that on the scale of this figure,

the spin orbit effects on the dispersion curve are rather modest compared to those in the

monolayer. For small wave vectors, with spin orbit coupling present, the dispersion curve of

the acoustic spin wave branch is fitted by the form Ω(Q‖) = 0.49 − 0.85Q‖ + 243Q 2 ‖(meV)

so at long wavelengths the influence of the Dzyaloshinskii Moriya interaction is more than

one order of magnitude smaller than it is in the monolayer.

If the magnetization is perpendicular to the surface, then symmetry considerations show

that there are no left/right asymmetries in the spin wave propagation characteristics. One

may see this as follows. Consider a wave vector ~Q‖ in the plane of the surface, which

also is perpendicular to the magnetization, and thus perpendicular to the long axis. The

reflection R in the plane perpendicular to the surface and which contains the magnetization

21

0 100 200 300 400 Ω (meV)

0

500

1000

1500

2000

2500

3000

Im χ +

,−

FIG. 8: (color online) For the bilayer and the case where the magnetization is perpendicular to the

surface (model B), and for Q‖ = 0.6A˚ −1 , we show spectral density function calculated for positive

values of Q‖ (continuous lines) and negative values of Q‖ (symbols). The black curve is the spectral

density for the inner layer, and the red curve is the outermost Fe layer.

simultaneously changes the sign of wave vector and the magnetization. If this is followed

by the time reversal operation T , then ~Q‖ remains reversed in sign but the magnetization

changes back to its original orientation. Thus the product RT leaves the system invariant

but transforms ~Q‖ into −~Q‖. The two propagation directions are then equivalent.

We illustrate this in Fig. 8 where, for Q‖ = 0.6A˚ −1 , where it is shown that the spectral

densities calculated for the two directions of propagation are identical, with spin orbit cou-

pling switched on. Model B, in which the magnetization is perpendicular to the surface, has

been used in these calculations. The spectral densities calculated for the two signs of Q‖

cannot be distinguished to within the numerical precision we use.

IV. CONCLUDING REMARKS

We have developed the formalism which allows one to include the influence of spin orbit

coupling on the spin excitations of ultrathin ferromagnets on semi infinite metallic substrates.

Our approach allows us to calculate the full dynamic susceptibility of the system, so as

illustrated by the calculations presented in section III we can examine the influence of spin

orbit coupling on the linewidth (or lifetime) of spin excitations, along with their oscillator

22

strength. As in previous work, we can then construct effective dispersion curves by following

peaks in the spectral density as a function of wave vector, without resort to calculations of

large numbers of very small distant neighbor exchange interactions. The results presented

in Fig. 4 are very similar to the experimental data reported in ref. 39, as discussed above,

though we see that in the bilayer the influence of the Dzyaloshinskii Moriya interaction is

considerably more modest than in the monolayer.

We will be exploring other issues in the near future. One interest in our minds is the

influence of spin orbit coupling on the spin pumping contributions to the ferromagnetic

resonance linewidth, as observed in ferromagnetic resonance (FMR) studies of ultrathin

films.43 It has been shown earlier44 that the methodology employed in the present paper

(without spin orbit coupling included) can be applied to the description of the spin pumping

contribution to the FMR linewidth, and in fact an excellent quantitative account of the

data on the Fe/Au(100) system was obtained. It is possible that for films grown on 4d

and 5d substrates that spin orbit coupling can influence the spin pumping relaxation rate

substantially. This will require calculations directed toward much thicker films than explored

here. The formalism we have developed and described in the present paper will allow such

studies in the future.

Acknowledgments

This research was supported by the U. S. Department of Energy, through grant No.

DE-FG03-84ER-45083. S. L. wishes to thank the Alexander von Humboldt Foundation

for a Feodor Lynen Fellowship. A.T.C. and R.B.M. acknowledge support from CNPq and

FAPERJ and A.B.K. was supported also by the CNPq, Brazil.

Appendix

In this Appendix we provide explicit expressions for the various quantities which enter

the equations displayed in Section II. While these expressions are unfortunately lengthy, it

will be useful for them to be given in full.

A (1) µν,µ′ν′(ll

′;mm′) = δl′mδνµ′〈c † lµ↑cm′ν′↑〉 − δlm′δµν′〈c

† mµ′↓cl′ν↓〉

23

A (2) µν,µ′ν′(ll

′;mm′) = −δlm′δµν′〈c † mµ↓cl′ν↑〉

A (3) µν,µ′ν′(ll

′;mm′) = δl′mδνµ′〈c † lµ↓cm′ν′↑〉

A (4) µν,µ′ν′(ll

′;mm′) = 0 (A.1)

The various expectation values in the equations above and those displayed below are

calculated from the single particle Greens functions once the self consistent ground state

parameters are determined. Then

B˜11µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul;µ′η,µν′〈c

† lη↓cl′ν↓〉δlmδlm′ − Ul′;µ′ν,ην′〈c

† lµ↑cl′η↑〉δl′mδl′m′

)

B˜12µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul′;µ′ν,ν′η〈c

† lµ↑cl′η↓〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↑cl′ν↓〉δlmδlm′+

+Ul;µ′η,µν′〈c † lη↑cl′ν↓〉δlmδlm′

)

B˜13µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul′;µ′ν,ν′η〈c

† lµ↑cl′η↓〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↑cl′ν↓〉δlmδlm′−

−Ul′;µ′ν,ην′〈c † lµ↑cl′η↓〉δl′mδl′m′

)

B˜14µν,µ′ν′(ll ′;mm′) = 0

(A.2)

B˜21µν,µ′ν′(ll ′;mm′) =

∑ η

Ul;µ′η,µν′〈c † lη↓cl′ν↑〉δlmδlm′

B˜22µν,µ′ν′(ll ′;mm′) =

∑ η

[ (Ul′;µ′ν,ν′η − Ul′;µ′ν,ην′) 〈c

† lµ↑cl′η↑〉δl′mδl′m′−

− (Ul;ηµ′,µν′ − Ul;µ′η,µν′) 〈c † lη↑cl′ν↑〉δlmδlm′

]

B˜23µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul;µ′ν,ν′η〈c

† lµ↑cl′η↑〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↑cl′ν↑〉δlmδl′m′

)

B˜24µν,µ′ν′(ll ′;mm′) = −

∑ η

Ul′;µ′ν,ην′〈c † lµ↑cl′η↓〉δl′mδl′m′

(A.3)

B˜31µν,µ′ν′(ll ′;mm′) = −

∑ η

Ul′;µ′ν,ην′〈c † lµ↓cl′η↑〉δl′mδl′m′

B˜32µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul;µ′ν,ν′η〈c

† lµ↓cl′η↓〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↓cl′ν↓〉δlmδl′m′

)

24

B˜33µν,µ′ν′(ll ′;mm′) =

∑ η

[ (Ul′;µ′ν,ν′η − Ul′;µ′ν,ην′) 〈c

† lµ↓cl′η↓〉δl′mδl′m′−

− (Ul;ηµ′,µν′ − Ul;µ′η,µν′) 〈c † lη↓cl′ν↓〉δlmδlm′

]

B˜34µν,µ′ν′(ll ′;mm′) =

∑ η

Ul;µ′η,µν′〈c † lη↑cl′ν↓〉δlmδlm′

(A.4)

B˜41µν,µ′ν′(ll ′;mm′) = 0

B˜42µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul′;µ′ν,ν′η〈c

† lµ↓cl′η↑〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↓cl′ν↑〉δlmδlm′−

−Ul′;µ′ν,ην′〈c † lµ↓cl′η↑〉δl′mδl′m′

)

B˜43µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul′;µ′ν,ν′η〈c

† lµ↓cl′η↑〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↓cl′ν↑〉δlmδlm′+

+Ul;µ′η,µν′〈c † lη↓cl′ν↑〉δlmδlm′

)

B˜44µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul;µ′η,µν′〈c

† lη↑cl′ν↑〉δlmδl′m − Ul′;µ′ν,ην′〈c

† lµ↓cl′η↓〉δl′mδl′m′

)

(A.5)

B11µν,µ′ν′(ll ′;mm′) = T˜ νν

′↓ l′m′ δlmδµµ′ − (T˜

µµ′↑ lm )

∗δl′m′δνν′

B12µν,µ′ν′(ll ′;mm′) = α∗l′;ν′νδlmδl′m′δµµ′

B13µν,µ′ν′(ll ′;mm′) = −α∗l;µµ′δlmδl′m′δνν′

B14µν,µ′ν′(ll ′;mm′) = 0 (A.6)

B21µν,µ′ν′(ll ′;mm′) = αl′;νν′δlmδl′m′δµµ′

B22µν,µ′ν′(ll ′;mm′) = T˜ νν

′↑ l′m′ δlmδµµ′ − (T˜

µµ′↑ lm )

∗δl′m′δνν′

B23µν,µ′ν′(ll ′;mm′) = 0

B24µν,µ′ν′(ll ′;mm′) = −α∗l;µµ′δlmδl′m′δνν′ (A.7)

B31µν,µ′ν′(ll ′;mm′) = −αl;µ′µδlmδl′m′δνν′

25

B32µν,µ′ν′(ll ′;mm′) = 0

B33µν,µ′ν′(ll ′;mm′) = T˜ νν

′↓ l′m′ δlmδµµ′ − (T˜

µµ′↓ lm )

∗δl′m′δνν′

B34µν,µ′ν′(ll ′;mm′) = α∗l′;ν′νδlmδl′m′δµµ′ (A.8)

B41µν,µ′ν′(ll ′;mm′) = 0

B42µν,µ′ν′(ll ′;mm′) = −αl;µ′µδlmδl′m′δνν′

B43µν,µ′ν′(ll ′;mm′) = αl′;νν′δlmδl′m′δµµ′

B11µν,µ′ν′(ll ′;mm′) = T˜ νν

′↑ l′m′ δlmδµµ′ − (T˜

µµ′↓ lm )

∗δl′m′δνν′ (A.9)

1 Rodrigo Arias and D. L. Mills, Phys. Rev. B 60, 7395 (1999); A. Azevedo, A. B. Oliveira, F. M.

de Aguiar, and S. M. Rezende, Phys. Rev. B 62, 5331 (2000); D. L. Mills and S. M. Rezende

in Spin Dynamics in Confined Structures II, (Springer Verlag, Heidelberg, 2002); P. Landeros,

Rodrigo E. Arias and D. L. Mills, Phys. Rev. B 77, 214405, (2008); J. Lindner, I. Barsukov, C.

Raeder, C. Hassel, O. Posth, R. Meckenstock, P. Landeros and D. L. Mills, Phys. Rev. B 80,

224421 (2009).

2 R. Urban, G. Woltersdorf and G. Heinrich, Phys. Rev. Lett. 87, 217204 (2001), Y. Tserkovnyak,

A. Brataas and G. E. W. Bauer, Phys. Rev. Lett. 88, 117601 (2002); Phys. Rev. B 66, 224403

(2002); D. L. Mills, Phys. Rev. B 68, 014419 (2003); M. Zwierzycki, Y. Tserkovnyak, P. J.

Kelly, A. Brataas, and G. E. W. Bauer, Phys. Rev. B 71, 064420 (2005); A. T. Costa, R. B.

Muniz and D. L. Mills, Phys. Rev. B 73, 054426 (2006).

3 S. M. Rezende, A. Azevedo, M. A. Lucena and F. M. Aguiar, Phys. Rev. B 63, 214416 (2001).

4 R. Vollmer, M. Etzkorn, P. S. Anil Kumar, H. Ibach and J. R. Kirschner, Phys. Rev. Lett.

91, 14720 (2003); for a theoretical discussion of the linewidth and dispersion observed in the

experiment, see A. T. Costa, R. B. Muniz Phys. Rev. B 70, 0544406 (2004).

5 Early papers are H. Tang, M. Plihal and D. L. Mills, Journal of Magnetism and Magnetic

Materials 187, 23 (1998); R. B. Muniz and D. L. Mills, Phys. Rev. B 66, 174417 (2002); A. T.

Costa, R. B. Muniz and D. L. Mills, Phys. Rev. B 68, 224435 (2003).

26

6 For an explicit comparison between the predictions of the Heisenberg model and that of an

itinerant electron description of the spin excitations, see Fig. 4 and Fig. 5 of the theoretical

paper cited in ref. 4.

7 S. Blgel, M. Weinert and P. H. Dederichs, Phys. Rev. Lett. 60, 1077 (1988).

8 S. Heinze, M. Bode, A. Kubetzka, O. Pietzsch, X. Nie, S. Blgel, and R. Wiesendanger, Science

311, 1805 (2000).

9 M. Bode, M. Heide, K. von Bergmann, P. Ferriani, S. Heinze, G. Bihlmayer, A. Kubetzka, O.

Pietzsch, S. Blgel and R. Wiesendanger, Nature 447, 190 (2007).

10 P. Lederer and A.Blandin, Phil. Mag. 14, 363 (1966).

11 See also A. T. Costa, R. B. Muniz and D. L. Mills, Phys. Rev. B 73, 214403 (2006).

12 L. Udvardi and L. Szunyogh, Phys. Rev. Lett. 102, 207204 (2009).

13 E. Abate and M. Asdente, Phys. Rev. 140, A1303 (1965).

14 A. T. Costa, Roberto Bechara Muniz and D. L. Mills, Phys. Rev. B 73, 054426 (2006).

15 Peter Fulde and Alan Luther, Phys. Rev. 175, 337 (1968).

16 See the discussion in the second paper cited in footnote 5.

17 See Chapter 6 and Chapter 7 and Appendix C of Group Theory and Quantum Mechanics, M.

Tinkham, (McGraw Hill, New York, 1964).

18 See D. M. Edwards, page 114 of Moment Formation in Solids, Vol. 117 of NATO Advanced

Study Institute, Series B: Physics (Plenum, New York, 1984), edited by W. J. L. Buyers.

19 See the second and third papers cited in footnote 5.

20 A. T. Costa, R. B. Muniz, J. X. Cao, R. Q. Wu and D. L. Mills, Phys. Rev. B 78, 054439

(2008).

21 J. Prokop, W. X. Tang, Y. Zhang, I Tudosa, T. R. F. Peixoto, Kh. Zakeri, and J. Kirschner,

Phys. Rev. Lett. 102, 177206 (2009).

22 W. X. Tang, Y. Zhang, I. Tudosa, J. Prokop, M. Etzkorn and J. Kirschner, Phys. Rev. Lett.

99, 087202 (2007).

23 See the discussion in A. T. Costa, R. B. Muniz and D. L. Mills, IEEE Transactions on Magnetics

44, 1974 (2008).

24 S. Krause, G. Herzog, T. Stapelfeldt, L. Berbil-Bautista, M. Bode E. Y. Vedmedenko, and R.

Wiesendanger, Phys. Rev. Lett. 103, 127202 (2009).

25 M. Pratzer, H. J. Elmers, M. Bode, O. Pietzsch, A. Kubetzka and R. Wiesendanger, Phys. Rev.

27

Lett. 88, 057201 (2002).

26 N. Papanikolaou, R. Zeller and P. H. Dederichs, J. Phys. Condens. Matter 14, 2799 (2002).

27 S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58, 1200 (1980).

28 A. I. Liechtenstein, M. I. Katsnelson, V. P. Antropov, and V. A. Gubanov, J. Magnetism and

Magnetic Materials 67, 65 (1987).

29 P. R. Peduto, S. Frota-Pessa and M. S. Methfessel, Phys. Rev. B 44, 13283 (1991).

30 S. Frota-Pessoˆa, Phys. Rev. B 46, 14570 (1992).

31 R. Haydock, Solid State Physics, edited by H. Ehrenreich, F. Seitz and D. Turnbull, (Academic

Press, New York, 1980), vol. 35.

32 O. K. Andersen, Phys. Rev. B 12, 3060 (1975).

33 O. K. Andersen, O. Jepsen and D. Gltzel, Highlights of Condensed Matter Theory, edited by

F. Bassani, F. Fumi, and M. P. Tosi, (North Holland, Amsterdam, 1985).

34 U. von Barth and L. A. Hedin, J. Phys. C. 5, 1629 (1972).

35 O. Eriksson, B. Johansson, R. C. Albers, A. M. Boring, and M. S. S. Brooks, Phys. Rev. B 42,

2707 (1990).

36 N. Beer and D. Pettifor, The Electronic Structure of Complex Systems, (Plenum Press, New

York, 1984).

37 M. Plihal and D. L. Mills, Phys. Rev. B 58, 14407 (1998).

38 M. Heide, G. Bihlmayer and S. Blgel, Phys. Rev. B 78, 140403 (2008).

39 Kh. Zakeri, Y. Zhang, J. Prokop, T. H. Chuang, N. Sakr, W. X. Tang and J. Kirschner, Phys.

Rev. Lett. 104, 137203 (2010).

40 B. M. Hall and D. L. Mills, Phys. Rev. B 34, 8318 (1986); Burl M. Hall, D. L. Mills, Mohammed

H. Mohammed and L. L. Kesmodel, Phys. Rev. B 38, 5856 (1988).

41 M. P. Gokhale, A. Ormeci and D. L. Mills, Phys. Rev. B 46, 8978 (1992).

42 B. M. Hall, S. Y. Tong and D. L. Mills, Phys. Rev. Lett. 50, 1277 (1983).

43 R. Urban, G. Woltersdorf, and B. Heinrich, Phys. Rev. Lett. 87, 217204 (2001).

44 A. T. Costa, R. B. Muniz and D. L. Mills, Phys. Rev. B 73, 054426 (2006).

45 A. Bergman, A. Taroni, L. Bergqvist, J. Hellsvik, B. Hrvarsson, and O. Eriksson, Phys. Rev. B

81, 144416 (2010).

28

iv :1

00 4.

30 66

v1 [

co nd

-m at.

me s-h

all ]

18 A

pr 20

10

Spin Orbit Coupling and Spin Waves in Ultrathin Ferromagnets:

The Spin Wave Rashba Effect

A. T. Costa,1 R. B. Muniz,1 S. Lounis,2 A. B. Klautau,3 and D. L. Mills2

1Instituto de F´ısica, Universidade Federal Fluminense, 24210-340 Nitero´i, RJ, Brasil.

2Department of Physics and Astronomy,

University of California Irvine, California, 92697, U. S. A.

3Departamento de Fisica, Universidade Federal do Para´, Bele´m, PA, Brazil.

Abstract

We present theoretical studies of the influence of spin orbit coupling on the spin wave excitations

of the Fe monolayer and bilayer on the W(110) surface. The Dzyaloshinskii-Moriya interaction is

active in such films, by virtue of the absence of reflection symmetry in the plane of the film. When

the magnetization is in plane, this leads to a linear term in the spin wave dispersion relation for

propagation across the magnetization. The dispersion relation thus assumes a form similar to

that of an energy band of an electron trapped on a semiconductor surfaces with Rashba coupling

active. We also show SPEELS response functions that illustrate the role of spin orbit coupling in

such measurements. In addition to the modifications of the dispersion relations for spin waves, the

presence of spin orbit coupling in the W substrate leads to a substantial increase in the linewidth

of the spin wave modes. The formalism we have developed applies to a wide range of systems,

and the particular system explored in the numerical calculations provides us with an illustration

of phenomena which will be present in other ultrathin ferromagnet/substrate combinations.

1

I. INTRODUCTION

The study of spin dynamics in ultrathin ferromagnets is of fundamental interest, since new

physics arises in these materials that has no counterpart in bulk magnetism. Examples are

provided by relaxation mechanisms evident in ferromagnetic resonance and Brillouin light

scattering studies,1–3 and also for the large wave vectors probed by spin polarized electron

loss spectroscopy (SPEELS).4 Of course, by now the remarkable impact of ultrathin film

structures on magnetic data storage is very well known, and other applications that exploit

spin dynamics in such materials are envisioned. Thus these issues are important from a

practical point of view as well as from that of fundamental physics.

Theoretical studies of the nature of spin waves in ultrathin films adsorbed on metal sub-

strates have been carried out for some years now, along with comparison with descriptions

provided with the Heisenberg model.5 In this paper, we extend the earlier theoretical treat-

ments to include the influence of spin orbit coupling on the spin wave spectrum of ultrathin

films. This extension is motivated by a most interesting discussion of the ground state of the

Mn monolayer on the W(110) surface. A nonrelativistic theoretical study of this system pre-

dicted that the ground state would be antiferromagnetic in character.7 This prediction was

confirmed by spin polarized scanning tunneling microscope studies of the system.8 However,

recent experimental STM data with a more sensitive instrument showed a more complex

ground state, wherein the ground state is in fact a spin density wave.9 One can construct

the new state by beginning with the antiferromagnet, and then superimposing on this a

long wavelength modulation on the direction of the moments on the lattice. The authors

of ref. 9 argued that the lack of reflection symmetry of the system in the plane of the film

activates the Dzyaloshinskii Moriya (DM) interaction, and the new state has its origin in this

interaction. They also presented relativistic and ab initio calculations that gave an excellent

account of the new data. The reflection symmetry is broken simply by the presence of the

substrate upon which the film is grown. This argument to us is most intriguing, since one

can then conclude that the DM interaction must be active in any ultrathin ferromagnet;

the substrate is surely always present. The DM interaction has its origin in the spin orbit

interaction, which of course is generally very weak in magnets that incorporate the 3d tran-

sition elements as the moment bearing entities. However, in the case of the Mn monolayer

on W(110) hybridization between the Mn 3d and the W 5d orbitals activates the very large

2

W spin orbit coupling, with the consequence that the strength of the DM interaction can

be substantial, as illustrated by the calculations presented in ref. 9. One may expect to see

substantial impact of the DM interaction in other ultrathin magnets grown on 5d substrates,

and possibly 4d substrates as well.

We have here another example of new physics present in ultrathin magnets that is not

encountered in the bulk form of the material from which the ultrathin structure is fabricated.

The purpose of this paper is to present our theoretical studies of spin orbit effects on spin

waves and also on the dynamic susceptibility of a much studied ultrathin film/substrate

combination, the Fe monolayer and bilayer on W(110). We find striking effects. For instance,

when the magnetization is in plane, as we shall see the DM interaction introduces a term

linear in wave vector in the dispersion relation of spin waves. Thus the uniform spin wave

mode at zero wave vector acquires a finite group velocity. We find this to be in the range

of 2 × 105 cm/sec for the Fe monolayer on W(110). Furthermore, left/right asymmetries

appear in the SPEELS response functions. Thus, we shall see that spin orbit coupling has

clear effects on the spin excitations of transition metal ultrathin ferromagnets grown on 5d

substrates.

We comment briefly on the philosophy of the approach used here, and in various earlier

publications.5 Numerous authors proceed as follows. One may generate a description of the

magnetic ground state of the adsorbed films by means of an electronic structure calculation

based on density functional theory. It is then possible to calculate, within the framework of

an adiabatic approximation, effective Heisenberg exchange integrals Jij between the magnetic

moments in unit cell i and unit cell j. These may be entered into a Heisenberg Hamilto-

nian, and then spin wave dispersion relations may be calculated through use of spin wave

theory. It has been known for decades10 that in the itinerant 3d magnets, effective exchange

interactions calculated in such a manner have very long range in real space. Thus, one must

include a very large number of distant neighbors in order to obtain converged results. This

is very demanding to do with high accuracy for the very numerous distant neighbors, since

the exchange interactions become very small as one moves out into distant neighbor shells.

At a more fundamental level, as noted briefly above, discussions in earlier publications

show that in systems such as we study here, the adiabatic approximation breaks down badly,

with qualitative consequences.5 First, spin wave modes of finite wave vector have very short

lifetimes, by virtue of decay into the continuum of particle hole pairs (Stoner excitations)

3

even at the absolute zero of temperature5,11 whereas in Heisenberg model descriptions their

lifetime is infinite. In multi layer films, the earlier calculations show that as a consequence of

the short lifetime, the spectrum of spin fluctuations at large wave vectors contains a single

broad feature which disperses with wave vector in a manner similar to that of a spin wave;

this is consistent with SPEELS data on an eight layer film of Co on Cu(100).4 This picture

stands in contrast to that offered by the Heisenberg model, in which a film of N layers has

N spin wave modes for each wave vector, and each mode has infinite lifetime.

The method developed earlier, and extended here to incorporate spin orbit coupling,

takes due account of the breakdown of the adiabatic approximation and also circumvents

the need to calculate effective exchange interactions in real space out to distant neighbor

shells. We work directly in wave vector space through study of the wave vector and frequency

dependent susceptibility discussed below, denoted as χ+,−( ~Q‖,Ω; l⊥, l ′ ⊥). The imaginary part

of this object, evaluated for l⊥ = l ′ ⊥ and considered as a function of frequency Ω for fixed

wave vector ~Q‖ provides us with the frequency spectrum of spin fluctuations on layer l⊥ for

the wave vector chosen. Spin waves appear as peaks in this function, very much as they do

in SPEELS data, and in a manner very similar to that used by experimentalist we extract

a dispersion relation for spin waves by following the wave vector dependence of the peak

frequency. We never need to resort to a real space summation procedure over large number

of neighbors, coupled by very tiny exchange couplings. The spin wave exchange stiffness

can be extracted either by fitting the small wave vector limit of the dispersion relation so

determined, or alternatively by utilizing an expression derived earlier5 which once again does

not require a summation in real space.

We comment on another feature of the present study. In earlier calculations,5,11,14 as in the

present paper, an empirical tight binding description forms the basis for our description of

the electronic structure. Within this approach, referred to as a multi band Hubbard model,

we can generate the wave vector and frequency dependent susceptibility for large systems.

In the earlier papers, effective tight binding parameters were extracted from bulk electronic

structure calculations. The present studies are based on tight binding parameters obtained

directly from a RS-LMTO-ASA calculation for the Fe/W(110) system. We also obtain tight

binding parameters by fitting KKR based electronic structure calculations for the ultrathin

film/substrate combinations of interest. We find that spin waves in the Fe/W(110) system

are quite sensitive to the empirical tight binding parameters which are employed, though as

4

we shall see the various descriptions provide very similar pictures of the one electron local

density of states.

We note that Udvardi and Szunyogh12 have also discussed the influence of spin orbit

coupling on the dispersion relation of spin waves in the Fe monolayer on W(110) within

the framework of the adiabatic approach discussed above, where exchange interactions and

other magnetic parameters are calculated in real space. We shall discuss a comparison with

our results and theirs below. There are differences. Most particularly, we note that in

Fig. 3, the authors of ref. 12 provide two dispersion curves for propagation perpendicular

to the magnetization, whereas in a film such as this with one spin per unit cell there can

be only one magnon branch. Additionally and very recently, Bergmann and coworkers45

investigated within an adiabatic approach finite temperature effects on the magnon spectrum

of Fe/W(110).

In section II, we comment on our means of introducing spin orbit coupling into the theory.

The results of our calculations are summarized in Section III and concluding remarks are

found in section IV.

II. CALCULATION OF THE DYNAMIC SUSCEPTIBILITY IN THE PRESENCE

OF SPIN ORBIT COUPLING

The formalism for including spin orbit coupling effects in our description of spin dynamics

is quite involved, so in this section we confine our attention to an outline of the key steps,

and an exposition of the overall structure of the theory. Our starting point is the multi band

Hubbard model of the system that was employed in our earlier study of spin dynamics in

ultrathin ferromagnets. The starting Hamiltonian is written as5

H = ∑ ij

∑ µνσ

T µν ij c

† iµσcjνσ +

1

2

∑ µνµ′ν′

∑ iσσ′

Ui;µν,µ′ν′c † iµσc

† iνσ′ciν′σ′ciµ′σ (1)

where i and j are site indices, σ, σ′ refer to spin, and µ, ν to the tight binding orbitals, nine in

number for each site, which are included in our treatment. The Coulomb interactions operate

only within the 3d orbitals on a given lattice site. The film, within which ferromagnetism

is driven by the Coulomb interactions, sits on a semi-infinite substrate within which the

Coulomb interaction is ignored.

5

In our empirical tight binding picture, the spin orbit interaction adds a term we write as

HSO = ∑ i

∑ µν

λi

2

[ Lzµν(c

† iµ↑ciν↑ − c

† iµciν↓) + L

+ µνc

† iµ↓ciν↑ + L

− µνc

† iµ↑ciν↓

] (2)

where ~L is the angular momentum operator, λi is the local spin-orbit coupling constant,

L± = Lx± iLy and Lαµν = 〈µ|L α|ν〉. We assume that the spin orbit interaction, present both

within the ferromagnetic film and the substrate, operates only within the 3d atomic orbitals.

A convenient tabulation of matrix elements of the orbital angular momentum operators is

found in ref. 13.

Information on the spin waves follows from the study of the spectral density of the trans-

verse dynamic susceptibility χ+,−( ~Q‖,Ω; l⊥, l ′ ⊥) as discussed above. From the text around

Eq. (1) of ref. 14, we see that this function describes the amplitude of the transverse spin

motion (the expectation value of the spin operator S+ in the layer labeled l⊥) to a fictitious

transverse magnetic field of frequency Ω and wave vector ~Q‖ parallel to the film surface that

is applied to layer l′⊥ of the sample. The spectral density, given by Im{χ+,−( ~Q‖,Ω; l⊥, l

′ ⊥)},

when multiplied by the Bose Einstein function n(Ω) = [exp(βΩ)−1]−1 is also the amplitude

of thermal spin fluctuations of wave vector ~Q‖ and frequency Ω in layer l⊥. We obtain in-

formation regarding the character (frequency, linewidth, and amplitude in layer l⊥) of spin

waves from the study of this function, as discussed earlier.5

Our previous analyses are based on the study of the dynamic susceptibility just described

through use of the random phase approximation (RPA) of many body theory. The Feynman

diagrams included in this method are the same as those incorporated into time dependent

density functional theory, though use of our Hubbard model allows us to solve the result-

ing equation easily once the very large array of irreducible particle hole propagators are

generated numerically.

Our task in the present paper is to extend the RPA treatment to incorporate spin orbit

coupling. The extension is non trivial. The quantity of interest, referred to in abbreviated

notation as χ+,−, may be expressed as a commutator of the spin operators S + and S− whose

precise definition is given earlier.5,12 With spin orbit coupling ignored, the RPA decoupling

procedure leads to a closed equation for χ+,−. When the RPA decoupling is carried out in

its presence, we are led to a sequence of four coupled equations which include new objects

we may refer to as χ−,−, χ↑,− and χ↓,−. The number of irreducible particle hole propagators

that must be computed likewise is increased by a factor of four. For a very simple version

6

of a one band Hubbard model, and for a very different purpose, Fulde and Luther carried

out an equivalent procedure many years ago15. In what follows, we provide a summary of

key steps along with expressions for the final set of equations.

To generate the equation of motion, we need the commutator of the operator S+µν(l, l ′) =

c † lµ↑cl′ν↓ with the Hamiltonian. One finds

[S+µν(l, l ′), HSO] =

1

2

∑ η

{λl′L + νηc

† lµ↑cl′η↑ − λlL

+ ηµc

† lη↓cl′ν↓ + λl′L

z νηc

† lµ↑cl′η↓ − λlL

z ηµc

† lη↑cl′ν↓}. (3)

The last two terms on the right hand side of Eq. 3 lead to terms in the equation of motion

which involve χ+,− whereas the first two terms couple us to the entities χ↑,− and χ↓,−.

When we write down the commutator of these new correlation functions with the spin orbit

Hamiltonian, we are led to terms which couple into the function χ−,− which is formed from

the commutator of two S− operators. In the absence of spin orbit coupling, a consequence of

spin rotation invariance of the Hamiltonian is that the three new functions just encountered

vanish. But they do not in its presence, and they must be incorporated into the analysis.

One then introduces the influence of the Coulomb interaction into the equation of motion,

and carries out an RPA decoupling of the resulting terms. The analysis is very lengthy, so

here we summarize only the structure that results from this procedure. Definitions of the

various quantities that enter are given in the Appendix. We express the equations of motion

in terms of a 4 × 4 matrix structure, where in schematic notation we let χ(1) = χ+,−,

χ(2) = χ↑,−, χ (3) = χ↓,− and χ

(4) = χ−,−. The four coupled equations then have the form

Ωχ(s) = A(s) + ∑ s′

(Bss ′

+ B˜ss ′

)χ(s ′) (4)

Each quantity in Eq. 4 has attached to it four orbital indices, and four site indices. To be

explicit, χ(2) = χ↑,− which enters Eq. (4) is formed from the commutator of the operator

c † lµ↑cl′ν↑ with c

† mµ′↓cm′ν′↑ and in full we denote this quantity as χ

(2) µν;µ′ν′(ll

′;mm′). The site

indices label the planes in the film, and we suppress reference to Ω and ~Q‖. The products

on the right hand side of Eq. 4 are matrix multiplications that involve these various indices.

For instance, the object ∑

s′ B ss′χ(s

′) is labeled by four orbital and four site indices so

[Bss ′

χ(s ′)]µν,µ′ν′(ll

′;mm′) = ∑ γδ

∑ nn′

Bss ′

µν,γδ(ll ′;nn′)χ

(s′) γδ,µ′ν′(nn

′;mm′). (5)

One proceeds by writing Eq. 4 in terms of the dynamic susceptibilities that characterize

the non-interacting system. These, referred to also as the irreducible particle hole propaga-

tors, are generated by evaluating the commutators which enter into the definition of χ(s) in

7

the non interacting ground state. These objects, denoted by χ(0s) obey a structure similar

to Eq. 4,

Ωχ(0s) = A(s) + ∑ s′

Bss ′

χ(s ′). (6)

It is then possible to relate χ(s) to χ(0s) through the relation, using four vector notation,

~χ(Ω) = ~χ(0)(Ω) + (Ω− B)−1B˜~χ(Ω). (7)

The matrix structure Γ ≡ (Ω−B)−1 may be generated from the definition of B, which may

be obtained from the equation of motion of the non-interacting susceptibility, Eq. 6. Then

B˜ follows from the equation of motion of the full susceptibility, as generated in the RPA.

One may solve Eq. 7

~χ(Ω) = [I − (Ω−B)−1B¯]−1~χ(0)(Ω), (8)

so our basic task is to compute the non interacting susceptibility matrix ~χ(0) and then carry

out the matrix inversion operation displayed in Eq. 8. For this we require the single particle

Greens functions (SPGFs) associated with our approach.

To generate the SPGFs, we set up an effective single particle Hamiltonian Hsp by intro-

ducing a mean field approximation for the Coulomb interaction. The general structure of

the single particle Hamiltonian is

Hsp = ∑ ij

∑ µνσ

T˜ µνσ ij c

† iµσcjνσ +

∑ i

∑ µν

{α∗i;µνc † iµ↓ciν↑ + αi;µνc

† iµ↑ciν↓} (9)

where the effective hopping integral T˜ µνσij contains the spin diagonal portion of the spin orbit

interaction, along with the mean field contributions from the Coulomb interaction. The form

we use for the latter is stated below. The coefficients in the spin flip terms are given by

αi;µν = λiL − µν −

∑ ηγ

Ui;ηµ,νγ〈c † iη↓ciγ↑〉. (10)

We then have the eigenvalue equation that generates the single particle eigenvalues and

eigenfunctions in the form Hsp|φs〉 = Es|φs〉; we can write this in the explicit form

∑ l

∑ ησ′

[ δσσ′ T˜

µησ′

il + δil(δσ′↓δσ↑α ∗ l;µη + δσ′↑δσ↓αl;µη)

] 〈lησ′|φs〉 = Es〈iµσ|φs〉. (11)

The single particle Greens function may be expressed in terms of the quantities that enter

Eq. 11. We have for this object the definition

Giµσ;jνσ′(t) = −iθ(t)〈{ciµσ(t), c † jνσ′(0)}〉 (12)

8

and one has the representation

Giµσ;jνσ′(Ω) = ∑ s

〈iµσ|φs〉〈φs|jνσ ′〉

Ω−Es + iη . (13)

These functions may be constructed directly from their equations of motion, which read,

after Fourier transforming with respect to time,

− ∑ l

∑ ησ′′

[ δσσ′′ T˜

µησ′′

il + δil(δσ′′↓δσ↑α ∗ l;µη + δσ′′↑δσ↓αl;µη)

] Glησ′′;jνσ′ + ΩGiµσ;jνσ′ = δσσ′δµνδij.

(14)

For the case where the substrate is semi infinite, our means of generating a numerical solution

to the hierarchy of equations stated in Eq. 14 has been discussed earlier. What remains is

to describe how the Coulomb interaction enters the effective hopping integrals T˜ µνσij that

appear in Eq. 9, Eq. 11 and Eq. 14.

There are, of course, a large number of Coulomb matrix elements in the original Hamil-

tonian, even if the Coulomb interactions are confined to within the 3d shell. Through the

use of group theory,17 the complete set of Coulomb matrix elements may be expressed in

terms of three parameters. These are given in Table I of the first cited paper in ref. 5.

In subsequent work, we have found that a much simpler structure18 nicely reproduces re-

sults obtained with the full three parameter form. We use the simpler one parameter form

here, for which Ui;µν,µ′ν′ = Uiδµν′δµ′ν . Then in the mean field approximation, the Coulomb

contribution to the single particle Hamiltonian assumes the form

H(C)sp = − ∑ i

Uimi

2

∑ µ

(c†iµ↑ciµ↑ − c † iµ↓ciµ↓) (15)

Here mi is the magnitude of the moment on site i. The Coulomb interactions Ui are non zero

only within the ultrathin ferromagnet, and the moments mi, determined self consistently,

vary from layer to layer when we consider multi layer ferromagnetic films.

It should be noted that when the Ansatz just described is employed in Eq. 10, the term

from the Coulomb interaction on the right hand side becomes proportional to the transverse

component of the moment located on site i and this vanishes identically. Thus, despite the

complexity introduced by the spin orbit coupling, when the simple one parameter Ansatz for

the Coulomb matrix elements is employed, one needs no parameters beyond the moment on

each layer in the self consistent loops that describe the ground state. In the present context,

this is an extraordinarily large savings in computational labor, and this will allow us to

9

address very large systems in the future. It is the case that certain off diagonal elements

such as 〈c†mµ′↓cl′ν↑〉 appear in the quantities defined in the Appendix. Notice, for example,

the expressions in Eqs. A.1. Once the ground state single particle Greens functions are

determined, such expectation values are readily computed.

III. RESULTS AND DISCUSSION

In earlier studies of Fe layers on W(110),19,20 as noted above, the electronic structure

was generated through use of tight binding parameters obtained from bulk electronic struc-

ture calculations. These calculations generate effective exchange interactions comparable

in magnitude to those found in the bulk transition metals,20 with the consequence that for

both monolayer Fe and bilayer Fe on W(110) the large wave vector spin waves generated

by theory are very much stiffer than found experimentally21,22 though it should be noted

that for the bilayer, the calculated value of the spin wave exchange stiffness is in excellent

accord with the data.23 Subsequent calculations which construct the spin wave dispersion

relation from adiabatic theory based on calculations of effective exchange integrals also gen-

erate spin waves for the monolayer substantially stiffer than found experimentally,12 though

they are softer than in our earlier work by a factor of two or so. We remark that it has

been suggested that the remarkably soft spin waves found experimentally may have origin

in carbon contamination of the monolayer and bilayer.20 We remark here that this can be

introduced during the SPEELS measurement. We note that the magnetic properties of Fe

monolayers grown on carbon free W(110)24 differ dramatically from those grown on surfaces

now known to be contaminated by carbon.25 In the former case, the domain walls have a

thickness of 2.15 nm,24 whereas in the latter circumstance very narrow walls with thickness

bounded from above by 0.6 nm are found.25 This suggests that the strength of the effective

exchange is very different in the two cases, with stiffer exchange in the carbon free samples.

The considerations of the previous paragraph have motivated us to carry out a series of

studies of the effective exchange in the Fe monolayer on W(110) within the framework of

three different electronic structure calculations. We find that although all three give local

density of states that are very similar, along with very similar energy bands when these are

examined, the intersite exchange interactions vary substantially. First, we have employed the

parameter set used earlier that is based on bulk electronic structures19,20 in new calculations

10

we call case A. In case B, we have employed an approach very similar to that used in ref. 12,

though in what follows our calculation of effective exchange integrals is non relativistic.

This is the Korringa Kohn Rostoker Greens Function (KKR-GF) method,26 which employs

the atomic sphere approximation and makes use of the Dyson equation G = g + gV G as

given in matrix notation. This allows us to calculate the Greens function G of an arbitrary

complex system given the perturbing potential V and the Greens function g of a reference

unperturbed system. Within the Local Spin Density Approximation (LSDA),27 We consider

a slab of five monolayers of W with the experimental lattice constant on top of which an Fe

monolayer is deposited and relaxed by -12.9%12 with respect to the W interlayer distance.

Angular momenta up to lmax = 3 were included in the Greens functions with a k mesh of

6400 points in the full two dimensional Brillouin zone. The effective exchange interactions

were calculated within the approximation of infinitesmal rotations28 that allows one to use

the magnetic force theorem. This states that the energy change due to infinitesmal rotations

in the moment directions can be calculated through the Kohn Sham eigenvalues.

Method C is the Real Space Linear-Muffin-Tin-Orbital approach as implemented, also,

in the atomic sphere approximation (RS-LMTO-ASA).29–33 Due to its linear scaling, this

method allows one to address the electronic structure of systems with a large number of

atoms for which the basic eigenvalue problem is solved in real space using the Haydock

recursion method. The Fe overlayer on the W(110) substrate was simulated by a large

bcc slab which contained ∼6800 atoms, arranged in 12 atomic planes parallel to the (110)

surface, with the experimental lattice parameter of bulk W. One empty sphere overlayer

is included, and self consistent potential parameters were obtained for the empty sphere

overlayer, the Fe monolayer, and the three W layers underneath using LSDA.34 For deeper

W layers we use bulk potential. Nine orbitals per site (the five 3d and 4 sp complex) were

used to describe the Fe valence band and the empty sphere overlayer, and for W the fully

occupied 4f orbitals were also included in the core. To evaluate the orbital moments we use

a scalar relativistic (SR) approach and include a spin orbit coupling term λ~L · ~S at each

variational step.35 In the recursion method the continued fraction has been terminated after

30 recursion levels with the Beer Pettifor terminator.36 The TB parameters so obtained are

inserted into our semi-empirical scheme and this allows us to generate the non interacting

susceptibilities which enter our full RPA description of the response of the structure.

In order to compare the electronic structures generated by the approaches just described,

11

we turn our attention to the local density of states for the majority and minority spins in

the adsorbed Fe monolayer. These are summarized in Fig. 1.

The local densities of states (LDOS) generated by the three sets of TB parameters have

approximately the same overall features, as we see from Fig. 1. The main differences appear

in the majority spin band, which overlaps the 5d states in the W substrate over a larger

energy range than the minority band. This is also true if we compare the LDOS generated

by the tight binding parameters extracted from the KKR electronic structure to the LDOS

obtained directly from the KKR calculations (red dashed line in Fig. 1b). The Fe-W hopping

parameters are indeed the least accurate portion of our parametrization scheme. In case A

we just used the Fe-Fe bulk parameters to describe the Fe-W hopping. In case B we extracted

TB parameters for Fe by fitting a KKR calculation of an unsupported Fe monolayer with a

lattice parameter matching that of the W substrate. For the Fe-W hopping we used the Fe

parameters obtained from the fitting, scaled to mimic the Fe-W distance relaxation. The

relaxation parameter was chosen to give the correct spin magnetic moment for the adsorbed

Fe monolayer. In case C all parameters were directly provided by the RS-LMTO-ASA code,

but in the DFT calculations the Fe-W distance was assumed to be equal to the distance

between W layers. Thus, the main difference between cases B and C is the treatment of the

mixing between Fe and W states and this is expected to affect more strongly the states that

occupy the same energy range.

As noted above, while the local density of states provided by the three approaches to

the electronic structure are quite similar as we see from Fig. 1 (and the same is true of the

electronic energy bands themselves if these are examined), the exchange interactions differ

substantially for the three cases. For the first, second and third neighbors we have (in meV)

42.5, 3.72 and 0.46 for model A, 28.7, -7.87 and 0.31, for model B and 11.23, -7.31 and 0.22

for model C. The authors of ref. 12 find 10.84, -3.34 and 3.64 for these exchange integrals.

We now turn to our studies of spin excitations in the Fe monolayer and Fe bilayer on

W(110) within the framework of the electronic structure generated through use of the ap-

proach in case C. We will discuss the influence of spin orbit coupling on both the trans-

verse wave vector dependent susceptibility though study of the spectral density function

A( ~Q‖,Ω; l⊥) = −Im{χ+,−( ~Q‖,Ω; l⊥)} discussed in section I. This function, for fixed wave

vector ~Q‖, when considered as a function of frequency Ω, describes the frequency spectrum

of the fluctuations of wave vector ~Q‖ of the transverse magnetic moment in layer l⊥ as noted

12

-6 -4 -2 0 2 E-EF (eV)

-2

0

2

4

LD O

S (st

ate s/e

V)

-2

0

2

4 LD

O S

(st ate

s/e V)

-2

0

2

4

LD O

S (st

ate s/e

V)

a)

b)

c)

FIG. 1: (color online) For the Fe monolayer adsorbed on W(110), we show the local density of

states in the Fe monolayer. The majority spin density of states is shown positive and the minority

spin density of states is negative. The zero of energy is at the Fermi energy. In (a), bulk electronic

structure parameters are used as in the second of the two papers cited in 5 (caseA). In (b), we

have the density of states generated by method B. The black curve is found by fitting the KKR

electronic energy band structure to tight binding parameters as described in the text, and the red

curve is calculated directly from the KKR calculation. In (c) we have the local density of states

generated by method C.

13

above. In the frequency regime where spin waves are encountered, this function is closely

related to (but not identical to) the response function probed in a SPEELS measurement.

In Fig. 2, for the Fe monolayer on W(110), we show the spectral density function cal-

culated for three values of | ~Q‖|, for propagation across the magnetization. Thus, the wave

vector is directed along the short axis in the surface. This is the direction probed in SPEELS

studies of the Fe monolayer on this surface.22 In each figure, we show three curves. The green

dashed curve is calculated with spin orbit coupling set to zero. We show only a single curve

for this case, because the spectral density is identical for the two directions of propagation

across the magnetization, + ~Q‖ and −~Q‖. When spin orbit coupling is switched on, for the

two directions just mentioned the response function is very different, as we see from the red

and black curve in the various panels. These spin wave frequencies, deduced from the peak in

the response functions as discussed in section I, differ for the two directions of propagation,

and also note that the peak intensities and linewidths differ as well. It is the absence of both

time reversal symmetry and reflection symmetry which renders + ~Q‖ and −~Q‖ inequivalent

for this direction of propagation. The system senses this breakdown of symmetry through

the spin orbit interaction. If one considers propagation parallel to the magnetization, the

asymmetries displayed in Fig. 2 are absent. The reason is that for this direction of propaga-

tion, reflection in the plane that is perpendicular to both the magnetization and the surface

is a good symmetry operation of the system, but takes + ~Q‖ into −~Q‖ thus rendering the two

directions equivalent. Recall, of course, that the magnetization is a pseudo vector in regard

to reflections. Notice how very broad the curves are for large wave vectors; the lifetime of

the spin waves is very short indeed.

As discussed in section I, we may construct a spin wave dispersion curve by plotting the

maxima in spectral density plots such as those illustrated in Fig. 2 as a function of wave

vector. We show dispersion relations constructed in this manner in Fig. 3, with spin orbit

coupling both present and absent. In Fig. 3a, and for propagation perpendicular to the

magnetization we show the dispersion curve so obtained for wave vectors throughout the

surface Brillouin zone, and in Fig. 3b we show its behavior for small wave vectors.

Let is first consider Fig, 3(a). Here the dispersion curve extends throughout the two

dimensional Brillouin zone. At the zone boundary, quite clearly the slope of the dispersion

curve does not vanish. In this direction of propagation, the nature of the point at the zone

boundary does not require the slope to vanish. What is most striking, clearly, is the anomaly

14

0 10 20 30 40 50 Ω (meV)

0

5×103

1×104

- Im

χ + ,−

(a)

0 50 100 150 200 Ω (meV)

0

500

1000

1500

- Im

χ + ,-

(b)

0 100 200 300 Ω (meV)

0

200

400

600

800

- Im

χ + ,-

(c)

FIG. 2: (color online) The spectral density functions A( ~Q‖,Ω; l⊥) evaluated in the Fe monolayer for

three values of the wave vector in the direction perpendicular to the magnetization. We have (a)

| ~Q‖| = 0.4A˚ −1 , (b) | ~Q‖| = 1.0A˚

−1 and (c) | ~Q‖| = 1.4A˚

−1 . The green curve (dashed) is calculated

with spin orbit coupling set to zero; the spectral density here is independent of the sign of ~Q‖. The

red and black curves are calculated with spin orbit coupling turned on. Now we see asymmetries

for propagation across the magnetization, with the red curve ~Q‖ directed from left to right and the

black curve from right to left.

in the vicinity of 1 A˚ −1 . This feature is evident in the calculation with spin orbit coupling

absent, and for positive values of the wave vector the feature becomes much more dramatic

when spin orbit coupling is switched on. Anomalies rather similar to those in the black

curve in Fig. 3(a) appear in the green dispersion curve found in Fig. 3 of ref. 12, though

these authors did not continue their calculation much beyond the 1 A˚ −1

regime. Our spin

waves are very much softer than theirs in this spectral region, notice. In Fig. 3 of ref. 12,

one finds two dispersion curves, one a mirror image of the second. Thus, these authors

15

-1.5 -1 -0.5 0 0.5 1 1.5

Q (A-1) 0

50

100

150

200

Ω

(m eV

) (a)

-0.4 -0.2 0 0.2 0.4

Q (A-1) 0

10

20

30

40

Ω

(m eV

)

(b)

FIG. 3: (color online) Spin wave dispersion relations constructed from peaks in the spectral den-

sity, for the Fe monolayer on W(110). The wave vector is in the direction perpendicular to the

magnetization. The red curve is constructed in the absence of spin orbit coupling, it is included in

the black curve.

display two spin wave frequencies for each wave vector. This surely is not correct. For a

structure with one atom per unit cell, there is one and only one spin wave mode for each

wave vector, though as discussed above for the structure explored here symmetry allow the

left/right asymmetry in the dispersion curve illustrated in our Fig. 3.

In Fig. 3b, again with spin orbit coupling switched on and off, we show an expanded

view of the dispersion curve for small wave vectors. With spin orbit interaction switched off,

at zero wave vector we see a zero frequency spin wave mode, as required by the Goldstone

theorem when the underlying Hamiltonian is form invariant under spin rotation. The curve

is also symmetrical, and is accurately fitted by the form Ω( ~Q‖) = 149Q 2 ‖ (meV), with the

wave vector in A˚ −1 , whereas with spin orbit coupling turned on the dispersion relation is

fitted by Ω( ~Q‖) = 3.4−11.8Q‖+143Q 2 ‖ (meV). Spin orbit coupling introduces an anisotropy

gap at Q‖ = 0, and most striking is the term linear in wave vector. This has its origin in the

Dzyaloshinskii Moriya interaction whose presence, as argued by the authors of ref. 9, has its

origin in the absence of both time reversal and inversion symmetry, for the adsorbed layer.

At long wavelengths, one may describe spin waves by classical long wavelength phe-

nomenology. The linear term in the dispersion curve has its origin in a term in the energy

density of the spin system of the form

VDM = −Γ ∫ dxdzSy(x, z)

∂Sx(x, z)

∂x (16)

16

Here Sα(x, z) is a spin density, the xz plane is parallel to the surface, and the magnetization

is parallel to the z direction.

One interesting feature of the spin wave mode whose dispersion relation is illustrated in

Fig. 3b is that at Q‖ = 0, the mode has a finite group velocity. The fit to the dispersion

curve gives this group velocity to be ∂Ω(~Q‖)

∂Q‖ ≈ 2×105cm/s , which is in the range of acoustic

phonon group velocities.

We turn now to our calculations of spin waves and the response functions for the Fe bilayer

on W(110). Let us first note that experimentally the orientation of the magnetization in

the bilayer appears to be dependent on the surface upon which the bilayer is grown. For

instance, when the bilayer is on the stepped W(110) surface, it is magnetized perpendicular

to the surface,25 a result in agreement with ab initio calculations of the anisotropy realized in

the epitaxial bilayer.38 However, in the SPEELS studies of spin excitations in the bilayer22,39

the magnetization is in plane. In our calculations, we find for model B the magnetization is

perpendicular to the surface, whereas in model C it lies in plane, along the long axis very

much as in the SPEELS experiments. The anisotropy in the bilayer is not particularly large,

on the order of 0.5 meV/Fe atom, and one sees from these results that it is a property quite

sensitive to the details of the electronic structure. The fact that model B and model C give

the two different stable orientation of the magnetization allows us to explore spin excitations

for the two different orientations of the magnetization.

We first turn our attention to the case where the magnetization lies in plane. The bilayer

has two spin wave modes, an acoustic mode for which the magnetization in the two planes

precesses in phase, and an optical mode for which they precess 180 degrees out of phase.

In Fig. 4, we show calculations of the dynamic susceptibility in the frequency range of the

acoustic mode for two values of the wave vector, Q‖ = +0.5A˚ − 1 and Q‖ = −0.5A˚

− 1. A

spin orbit induced left right asymmetry is clearly evident both in the peak frequency and

the height of the feature. Very recently, beautiful measurements of spin orbit asymmetries

in the Fe bilayer have appeared,39 and the results of our Fig. 4 are to be compared with

Fig. 3 of ref. 39. Theory and experiment are very similar, both in regard to the intensity

asymmetry and also the spin orbit induced frequency shift, though our calculated spin wave

frequencies are a little stiffer than those found experimentally.

As remarked above, in Fig. 4 we show only the acoustical spin wave mode frequency

regime. In Fig. 5, for the spectral densities in the innermost layer (upper panel) and the

17

0 50 100 Ω (meV)

0

500

1000

1500

2000

2500

3000

- Im

χ + ,−

FIG. 4: (color online) The spectral density in the innermost layer, in the acoustic spin wave regime,

for wave vectors of Q‖ = +0.5A˚ −1

(black curve) and Q‖ = −0.5A˚ −1

(red curve). Model C has

been used for the calculation. In the ground state, the magnetization lies in plane along the long

axis.

outermost layer (lower panel) we show the spectral densities for the entire spin wave regime,

including the region where the optical spin wave is found. It is clear that the spin orbit

induced frequency shifts are largest for the optical mode which, unfortunately is not observed

in the experiments.39

In Fig. 6 for a sequence of wave vectors, all chosen positive, we show a sequence of spectra

calculated for the entire frequency range so both the acoustic and optical spin wave feature

are displayed. The black curves show the spectral density of the innermost Fe layer, and

the red curves are for the outer layer. The optical spin wave mode, not evident in the data,

shows clearly in these figures. Notice that for wave vectors greater than 1A˚ −1

the acoustical

mode is localized in the outer layer and the optical mode is localized on the inner layer.

The optical mode is very much broader than the acoustical mode at large wave vectors, by

virtue of the strong coupling to the electron hole pairs in the W 5d bands.

An interesting issue is the absence of the optical mode from the SPEELS spectra reported

18

-1500

-1000

-500

0

- Im

χ + ,−

0 50 100 150 200 250 Ω (meV)

-2000

-1000

0

- Im

χ + ,−

FIG. 5: (color online) For the wave vector Q‖ = 0.5A˚ −1

we show the spectral densities in the

innermost Fe layer (upper panel) and in the outermost layer (lower panel) for the Fe bilayer on

W(110). The figure includes the optical spin wave feature. As in Fig. 4, the black curve is calculated

for Q‖ positive, and the red curves are for Q‖ negative. The calculations employ model C.

in refs. 22 and 39. We note that these spectra are taken with only two beam energies, 4 eV

and 6.75 eV. At such very low energies, the beam electron will sample both Fe layers,

so the SPEELS signal will be a coherent superposition of electron waves backscattered

from each layer; the excitation process involves coherent excitation of both layers by the

incident electron. As a consequence of the 180 degree phase difference in spin motions

associated with the two modes it is quite possible, indeed even probable, that for energies

where the acoustical mode is strong the intensity of the optical mode is weak, by virtue

of quantum interference effects in the excitation scattering amplitude. In earlier studies

of surface phonons, it is well documented that on surfaces where two surface phonons of

different polarization exist for the same wave vector, one can be silent and one active in

electron loss spectroscopy.40 It would require a full multiple scattering analysis of the spin

wave excitation process to explore this theoretically. While earlier41 calculations that address

SPEELS excitation of spin waves described by the Heisenberg model could be adapted for

19

0 2000 4000 6000

Im χ +

,−

0

1000

2000 Im

χ + ,−

0

500

1000

Im χ +

,−

0

500

1000

1500

Im χ +

,−

0 500

1000 1500

Im χ +

,−

0 100 200 300 400 500 Ω (meV)

0

500

1000

Im χ +

,−

Qy=0.4 A -1

Qy=0.6 A -1

Qy=0.8 A -1

Qy=1.0 A -1

Qy=1.2 A -1

Qy=1.4 A -1

FIG. 6: (color online) For the Fe bilayer and for several values of the wave vector (all positive), we

show the spectral density functions for the innermost layer adjacent to the substrate (black curve)

and those for the outer layer of the film. The calculations employ model C.

20

-1.5 -1 -0.5 0 0.5 1 1.5 Q (A-1)

0

100

200

300

Ω

(m eV

)

FIG. 7: (color online) For the Fe bilayer with magnetization in plane, we show the spin wave

dispersion curves calculated with spin orbit coupling (black points) and without spin orbit coupling

(red curves). Model C has been employed for these calculations.

this purpose, in principle, a problem is that at such low beam energies it is necessary to take

due account of image potential effects to obtain meaningful results.42 This is very difficult

to do without considerable information on the electron reflectivity of the surface.42 It would

be of great interest to see experimental SPEELS studies of the Fe bilayer with a wider range

of beam energies to search for the optical mode, if this were possible.

In Fig. 7, we show dispersion curves for the optical and acoustic spin wave branches for

the bilayer. The magnetization lies in plane, and one can see that on the scale of this figure,

the spin orbit effects on the dispersion curve are rather modest compared to those in the

monolayer. For small wave vectors, with spin orbit coupling present, the dispersion curve of

the acoustic spin wave branch is fitted by the form Ω(Q‖) = 0.49 − 0.85Q‖ + 243Q 2 ‖(meV)

so at long wavelengths the influence of the Dzyaloshinskii Moriya interaction is more than

one order of magnitude smaller than it is in the monolayer.

If the magnetization is perpendicular to the surface, then symmetry considerations show

that there are no left/right asymmetries in the spin wave propagation characteristics. One

may see this as follows. Consider a wave vector ~Q‖ in the plane of the surface, which

also is perpendicular to the magnetization, and thus perpendicular to the long axis. The

reflection R in the plane perpendicular to the surface and which contains the magnetization

21

0 100 200 300 400 Ω (meV)

0

500

1000

1500

2000

2500

3000

Im χ +

,−

FIG. 8: (color online) For the bilayer and the case where the magnetization is perpendicular to the

surface (model B), and for Q‖ = 0.6A˚ −1 , we show spectral density function calculated for positive

values of Q‖ (continuous lines) and negative values of Q‖ (symbols). The black curve is the spectral

density for the inner layer, and the red curve is the outermost Fe layer.

simultaneously changes the sign of wave vector and the magnetization. If this is followed

by the time reversal operation T , then ~Q‖ remains reversed in sign but the magnetization

changes back to its original orientation. Thus the product RT leaves the system invariant

but transforms ~Q‖ into −~Q‖. The two propagation directions are then equivalent.

We illustrate this in Fig. 8 where, for Q‖ = 0.6A˚ −1 , where it is shown that the spectral

densities calculated for the two directions of propagation are identical, with spin orbit cou-

pling switched on. Model B, in which the magnetization is perpendicular to the surface, has

been used in these calculations. The spectral densities calculated for the two signs of Q‖

cannot be distinguished to within the numerical precision we use.

IV. CONCLUDING REMARKS

We have developed the formalism which allows one to include the influence of spin orbit

coupling on the spin excitations of ultrathin ferromagnets on semi infinite metallic substrates.

Our approach allows us to calculate the full dynamic susceptibility of the system, so as

illustrated by the calculations presented in section III we can examine the influence of spin

orbit coupling on the linewidth (or lifetime) of spin excitations, along with their oscillator

22

strength. As in previous work, we can then construct effective dispersion curves by following

peaks in the spectral density as a function of wave vector, without resort to calculations of

large numbers of very small distant neighbor exchange interactions. The results presented

in Fig. 4 are very similar to the experimental data reported in ref. 39, as discussed above,

though we see that in the bilayer the influence of the Dzyaloshinskii Moriya interaction is

considerably more modest than in the monolayer.

We will be exploring other issues in the near future. One interest in our minds is the

influence of spin orbit coupling on the spin pumping contributions to the ferromagnetic

resonance linewidth, as observed in ferromagnetic resonance (FMR) studies of ultrathin

films.43 It has been shown earlier44 that the methodology employed in the present paper

(without spin orbit coupling included) can be applied to the description of the spin pumping

contribution to the FMR linewidth, and in fact an excellent quantitative account of the

data on the Fe/Au(100) system was obtained. It is possible that for films grown on 4d

and 5d substrates that spin orbit coupling can influence the spin pumping relaxation rate

substantially. This will require calculations directed toward much thicker films than explored

here. The formalism we have developed and described in the present paper will allow such

studies in the future.

Acknowledgments

This research was supported by the U. S. Department of Energy, through grant No.

DE-FG03-84ER-45083. S. L. wishes to thank the Alexander von Humboldt Foundation

for a Feodor Lynen Fellowship. A.T.C. and R.B.M. acknowledge support from CNPq and

FAPERJ and A.B.K. was supported also by the CNPq, Brazil.

Appendix

In this Appendix we provide explicit expressions for the various quantities which enter

the equations displayed in Section II. While these expressions are unfortunately lengthy, it

will be useful for them to be given in full.

A (1) µν,µ′ν′(ll

′;mm′) = δl′mδνµ′〈c † lµ↑cm′ν′↑〉 − δlm′δµν′〈c

† mµ′↓cl′ν↓〉

23

A (2) µν,µ′ν′(ll

′;mm′) = −δlm′δµν′〈c † mµ↓cl′ν↑〉

A (3) µν,µ′ν′(ll

′;mm′) = δl′mδνµ′〈c † lµ↓cm′ν′↑〉

A (4) µν,µ′ν′(ll

′;mm′) = 0 (A.1)

The various expectation values in the equations above and those displayed below are

calculated from the single particle Greens functions once the self consistent ground state

parameters are determined. Then

B˜11µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul;µ′η,µν′〈c

† lη↓cl′ν↓〉δlmδlm′ − Ul′;µ′ν,ην′〈c

† lµ↑cl′η↑〉δl′mδl′m′

)

B˜12µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul′;µ′ν,ν′η〈c

† lµ↑cl′η↓〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↑cl′ν↓〉δlmδlm′+

+Ul;µ′η,µν′〈c † lη↑cl′ν↓〉δlmδlm′

)

B˜13µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul′;µ′ν,ν′η〈c

† lµ↑cl′η↓〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↑cl′ν↓〉δlmδlm′−

−Ul′;µ′ν,ην′〈c † lµ↑cl′η↓〉δl′mδl′m′

)

B˜14µν,µ′ν′(ll ′;mm′) = 0

(A.2)

B˜21µν,µ′ν′(ll ′;mm′) =

∑ η

Ul;µ′η,µν′〈c † lη↓cl′ν↑〉δlmδlm′

B˜22µν,µ′ν′(ll ′;mm′) =

∑ η

[ (Ul′;µ′ν,ν′η − Ul′;µ′ν,ην′) 〈c

† lµ↑cl′η↑〉δl′mδl′m′−

− (Ul;ηµ′,µν′ − Ul;µ′η,µν′) 〈c † lη↑cl′ν↑〉δlmδlm′

]

B˜23µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul;µ′ν,ν′η〈c

† lµ↑cl′η↑〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↑cl′ν↑〉δlmδl′m′

)

B˜24µν,µ′ν′(ll ′;mm′) = −

∑ η

Ul′;µ′ν,ην′〈c † lµ↑cl′η↓〉δl′mδl′m′

(A.3)

B˜31µν,µ′ν′(ll ′;mm′) = −

∑ η

Ul′;µ′ν,ην′〈c † lµ↓cl′η↑〉δl′mδl′m′

B˜32µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul;µ′ν,ν′η〈c

† lµ↓cl′η↓〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↓cl′ν↓〉δlmδl′m′

)

24

B˜33µν,µ′ν′(ll ′;mm′) =

∑ η

[ (Ul′;µ′ν,ν′η − Ul′;µ′ν,ην′) 〈c

† lµ↓cl′η↓〉δl′mδl′m′−

− (Ul;ηµ′,µν′ − Ul;µ′η,µν′) 〈c † lη↓cl′ν↓〉δlmδlm′

]

B˜34µν,µ′ν′(ll ′;mm′) =

∑ η

Ul;µ′η,µν′〈c † lη↑cl′ν↓〉δlmδlm′

(A.4)

B˜41µν,µ′ν′(ll ′;mm′) = 0

B˜42µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul′;µ′ν,ν′η〈c

† lµ↓cl′η↑〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↓cl′ν↑〉δlmδlm′−

−Ul′;µ′ν,ην′〈c † lµ↓cl′η↑〉δl′mδl′m′

)

B˜43µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul′;µ′ν,ν′η〈c

† lµ↓cl′η↑〉δl′mδl′m′ − Ul;ηµ′,µν′〈c

† lη↓cl′ν↑〉δlmδlm′+

+Ul;µ′η,µν′〈c † lη↓cl′ν↑〉δlmδlm′

)

B˜44µν,µ′ν′(ll ′;mm′) =

∑ η

( Ul;µ′η,µν′〈c

† lη↑cl′ν↑〉δlmδl′m − Ul′;µ′ν,ην′〈c

† lµ↓cl′η↓〉δl′mδl′m′

)

(A.5)

B11µν,µ′ν′(ll ′;mm′) = T˜ νν

′↓ l′m′ δlmδµµ′ − (T˜

µµ′↑ lm )

∗δl′m′δνν′

B12µν,µ′ν′(ll ′;mm′) = α∗l′;ν′νδlmδl′m′δµµ′

B13µν,µ′ν′(ll ′;mm′) = −α∗l;µµ′δlmδl′m′δνν′

B14µν,µ′ν′(ll ′;mm′) = 0 (A.6)

B21µν,µ′ν′(ll ′;mm′) = αl′;νν′δlmδl′m′δµµ′

B22µν,µ′ν′(ll ′;mm′) = T˜ νν

′↑ l′m′ δlmδµµ′ − (T˜

µµ′↑ lm )

∗δl′m′δνν′

B23µν,µ′ν′(ll ′;mm′) = 0

B24µν,µ′ν′(ll ′;mm′) = −α∗l;µµ′δlmδl′m′δνν′ (A.7)

B31µν,µ′ν′(ll ′;mm′) = −αl;µ′µδlmδl′m′δνν′

25

B32µν,µ′ν′(ll ′;mm′) = 0

B33µν,µ′ν′(ll ′;mm′) = T˜ νν

′↓ l′m′ δlmδµµ′ − (T˜

µµ′↓ lm )

∗δl′m′δνν′

B34µν,µ′ν′(ll ′;mm′) = α∗l′;ν′νδlmδl′m′δµµ′ (A.8)

B41µν,µ′ν′(ll ′;mm′) = 0

B42µν,µ′ν′(ll ′;mm′) = −αl;µ′µδlmδl′m′δνν′

B43µν,µ′ν′(ll ′;mm′) = αl′;νν′δlmδl′m′δµµ′

B11µν,µ′ν′(ll ′;mm′) = T˜ νν

′↑ l′m′ δlmδµµ′ − (T˜

µµ′↓ lm )

∗δl′m′δνν′ (A.9)

1 Rodrigo Arias and D. L. Mills, Phys. Rev. B 60, 7395 (1999); A. Azevedo, A. B. Oliveira, F. M.

de Aguiar, and S. M. Rezende, Phys. Rev. B 62, 5331 (2000); D. L. Mills and S. M. Rezende

in Spin Dynamics in Confined Structures II, (Springer Verlag, Heidelberg, 2002); P. Landeros,

Rodrigo E. Arias and D. L. Mills, Phys. Rev. B 77, 214405, (2008); J. Lindner, I. Barsukov, C.

Raeder, C. Hassel, O. Posth, R. Meckenstock, P. Landeros and D. L. Mills, Phys. Rev. B 80,

224421 (2009).

2 R. Urban, G. Woltersdorf and G. Heinrich, Phys. Rev. Lett. 87, 217204 (2001), Y. Tserkovnyak,

A. Brataas and G. E. W. Bauer, Phys. Rev. Lett. 88, 117601 (2002); Phys. Rev. B 66, 224403

(2002); D. L. Mills, Phys. Rev. B 68, 014419 (2003); M. Zwierzycki, Y. Tserkovnyak, P. J.

Kelly, A. Brataas, and G. E. W. Bauer, Phys. Rev. B 71, 064420 (2005); A. T. Costa, R. B.

Muniz and D. L. Mills, Phys. Rev. B 73, 054426 (2006).

3 S. M. Rezende, A. Azevedo, M. A. Lucena and F. M. Aguiar, Phys. Rev. B 63, 214416 (2001).

4 R. Vollmer, M. Etzkorn, P. S. Anil Kumar, H. Ibach and J. R. Kirschner, Phys. Rev. Lett.

91, 14720 (2003); for a theoretical discussion of the linewidth and dispersion observed in the

experiment, see A. T. Costa, R. B. Muniz Phys. Rev. B 70, 0544406 (2004).

5 Early papers are H. Tang, M. Plihal and D. L. Mills, Journal of Magnetism and Magnetic

Materials 187, 23 (1998); R. B. Muniz and D. L. Mills, Phys. Rev. B 66, 174417 (2002); A. T.

Costa, R. B. Muniz and D. L. Mills, Phys. Rev. B 68, 224435 (2003).

26

6 For an explicit comparison between the predictions of the Heisenberg model and that of an

itinerant electron description of the spin excitations, see Fig. 4 and Fig. 5 of the theoretical

paper cited in ref. 4.

7 S. Blgel, M. Weinert and P. H. Dederichs, Phys. Rev. Lett. 60, 1077 (1988).

8 S. Heinze, M. Bode, A. Kubetzka, O. Pietzsch, X. Nie, S. Blgel, and R. Wiesendanger, Science

311, 1805 (2000).

9 M. Bode, M. Heide, K. von Bergmann, P. Ferriani, S. Heinze, G. Bihlmayer, A. Kubetzka, O.

Pietzsch, S. Blgel and R. Wiesendanger, Nature 447, 190 (2007).

10 P. Lederer and A.Blandin, Phil. Mag. 14, 363 (1966).

11 See also A. T. Costa, R. B. Muniz and D. L. Mills, Phys. Rev. B 73, 214403 (2006).

12 L. Udvardi and L. Szunyogh, Phys. Rev. Lett. 102, 207204 (2009).

13 E. Abate and M. Asdente, Phys. Rev. 140, A1303 (1965).

14 A. T. Costa, Roberto Bechara Muniz and D. L. Mills, Phys. Rev. B 73, 054426 (2006).

15 Peter Fulde and Alan Luther, Phys. Rev. 175, 337 (1968).

16 See the discussion in the second paper cited in footnote 5.

17 See Chapter 6 and Chapter 7 and Appendix C of Group Theory and Quantum Mechanics, M.

Tinkham, (McGraw Hill, New York, 1964).

18 See D. M. Edwards, page 114 of Moment Formation in Solids, Vol. 117 of NATO Advanced

Study Institute, Series B: Physics (Plenum, New York, 1984), edited by W. J. L. Buyers.

19 See the second and third papers cited in footnote 5.

20 A. T. Costa, R. B. Muniz, J. X. Cao, R. Q. Wu and D. L. Mills, Phys. Rev. B 78, 054439

(2008).

21 J. Prokop, W. X. Tang, Y. Zhang, I Tudosa, T. R. F. Peixoto, Kh. Zakeri, and J. Kirschner,

Phys. Rev. Lett. 102, 177206 (2009).

22 W. X. Tang, Y. Zhang, I. Tudosa, J. Prokop, M. Etzkorn and J. Kirschner, Phys. Rev. Lett.

99, 087202 (2007).

23 See the discussion in A. T. Costa, R. B. Muniz and D. L. Mills, IEEE Transactions on Magnetics

44, 1974 (2008).

24 S. Krause, G. Herzog, T. Stapelfeldt, L. Berbil-Bautista, M. Bode E. Y. Vedmedenko, and R.

Wiesendanger, Phys. Rev. Lett. 103, 127202 (2009).

25 M. Pratzer, H. J. Elmers, M. Bode, O. Pietzsch, A. Kubetzka and R. Wiesendanger, Phys. Rev.

27

Lett. 88, 057201 (2002).

26 N. Papanikolaou, R. Zeller and P. H. Dederichs, J. Phys. Condens. Matter 14, 2799 (2002).

27 S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58, 1200 (1980).

28 A. I. Liechtenstein, M. I. Katsnelson, V. P. Antropov, and V. A. Gubanov, J. Magnetism and

Magnetic Materials 67, 65 (1987).

29 P. R. Peduto, S. Frota-Pessa and M. S. Methfessel, Phys. Rev. B 44, 13283 (1991).

30 S. Frota-Pessoˆa, Phys. Rev. B 46, 14570 (1992).

31 R. Haydock, Solid State Physics, edited by H. Ehrenreich, F. Seitz and D. Turnbull, (Academic

Press, New York, 1980), vol. 35.

32 O. K. Andersen, Phys. Rev. B 12, 3060 (1975).

33 O. K. Andersen, O. Jepsen and D. Gltzel, Highlights of Condensed Matter Theory, edited by

F. Bassani, F. Fumi, and M. P. Tosi, (North Holland, Amsterdam, 1985).

34 U. von Barth and L. A. Hedin, J. Phys. C. 5, 1629 (1972).

35 O. Eriksson, B. Johansson, R. C. Albers, A. M. Boring, and M. S. S. Brooks, Phys. Rev. B 42,

2707 (1990).

36 N. Beer and D. Pettifor, The Electronic Structure of Complex Systems, (Plenum Press, New

York, 1984).

37 M. Plihal and D. L. Mills, Phys. Rev. B 58, 14407 (1998).

38 M. Heide, G. Bihlmayer and S. Blgel, Phys. Rev. B 78, 140403 (2008).

39 Kh. Zakeri, Y. Zhang, J. Prokop, T. H. Chuang, N. Sakr, W. X. Tang and J. Kirschner, Phys.

Rev. Lett. 104, 137203 (2010).

40 B. M. Hall and D. L. Mills, Phys. Rev. B 34, 8318 (1986); Burl M. Hall, D. L. Mills, Mohammed

H. Mohammed and L. L. Kesmodel, Phys. Rev. B 38, 5856 (1988).

41 M. P. Gokhale, A. Ormeci and D. L. Mills, Phys. Rev. B 46, 8978 (1992).

42 B. M. Hall, S. Y. Tong and D. L. Mills, Phys. Rev. Lett. 50, 1277 (1983).

43 R. Urban, G. Woltersdorf, and B. Heinrich, Phys. Rev. Lett. 87, 217204 (2001).

44 A. T. Costa, R. B. Muniz and D. L. Mills, Phys. Rev. B 73, 054426 (2006).

45 A. Bergman, A. Taroni, L. Bergqvist, J. Hellsvik, B. Hrvarsson, and O. Eriksson, Phys. Rev. B

81, 144416 (2010).

28

Comments