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