First principles theory of chiral dichroism in electron microscopy applied to 3d
ferromagnets
Ján Rusz∗
Institute of Physics, Academy of Sciences of the Czech Republic,
Cukrovarnická 10, 162 53 Prague 6, Czech Republic†
Stefano Rubino and Peter Schattschneider
Institute for Solid State Physics, Vienna University of Technology,
Wiedner Hauptstrasse 8-10/138, A-1040 Vienna, Austria
(Dated: August 2, 2021)
Recently it was demonstrated (Schattschneider et al., Nature 441 (2006), 486), that an analogue of
the X-ray magnetic circular dichroism (XMCD) experiment can be performed with the transmission
electron microscope (TEM). The new phenomenon has been named energy-loss magnetic chiral
dichroism (EMCD). In this work we present a detailed ab initio study of the chiral dichroism in the
Fe, Co and Ni transition elements. We discuss the methods used for the simulations together with
the validity and accuracy of the treatment, which can, in principle, apply to any given crystalline
specimen. The dependence of the dichroic signal on the sample thickness, accuracy of the detector
position and the size of convergence and collection angles is calculated.
PACS numbers:
Keywords: density functional theory, chiral dichroism, transmission electron microscopy, dynamical diffrac-
tion theory
I. INTRODUCTION
The analogy between X-ray absorption spectroscopy
(XAS) and electron energy loss spectroscopy (EELS) has
been recognized long ago1,2. The role of the polarization
vector ε in XAS is similar to the role of the wave vector
transfer q in EELS. This has made feasible the detection
of linear dichroism in the TEM. However the counter-
part of X-ray magnetic circular dichroism (XMCD)3,4,5
experiments with electron probes was thought to be tech-
nically impossible due to the low intensity of existing
spin polarized electron sources. XMCD is an important
technique providing atom-specific information about the
magnetic properties of materials. Particularly the near
edge spectra, where a well localized strongly bound elec-
tron with l 6= 0 is excited to an unoccupied band state,
allow to measure spin and orbital moments. Soon after
the proposal of an experimental setup for detection of cir-
cular dichroism using a standard non-polarized electron
beam in the TEM6 it was demonstrated that such ex-
periments (called energy-loss magnetic chiral dichroism,
EMCD) are indeed possible7. This novel technique is of
considerable interest for nanomagnetism and spintron-
ics according to the high spatial resolution of the TEM.
However, its optimization involves many open questions.
In this work we provide theoretical ab initio predictions
of the dependence of the dichroic signal in the EMCD
experiment on several experimental conditions, such as
sample thickness, detector placement and the finite size
of convergence and collection angles. This information
should help to optimize the experimental geometry in
order to maximize the signal to noise ratio.
The structure of this work is as follows: In section II
we first describe the computational approach based on
the dynamical diffraction theory and electronic structure
calculations. We also discuss the validity of several ap-
proximations for the mixed dynamic form factor. In sec-
tion III we study the dependence of the dichroic signal of
bcc-Fe, hcp-Co and fcc-Ni on various experimental con-
ditions. This section is followed by a concluding section
summarizing the most important findings.
II. METHOD OF CALCULATION
We will follow the derivations of the double differential
scattering cross-section (DDSCS) presented in Refs. 8,9.
Within the first-order Born approximation10 DDSCS is
written as
S(q, E)
S(q, E) =
|〈i|eiq·R̂|f〉|2δ(Ef − Ei − E) (2)
where q = kf−k0 is the difference (wave vector transfer)
between final wave vector kf and initial wave vector k0
of the fast electron; γ = 1/
1− v2/c2 is a relativistic
factor and a0 is Bohr radius. The S(q, E) is the so called
dynamic form factor (DFF)11.
This equation is valid only if the initial and final wave
functions of the fast electron are plane waves. In the
crystal the full translation symmetry is broken and as
a result, the electron wave function becomes a superpo-
sition of Bloch waves, which reflects the discrete trans-
lation symmetry. Each Bloch wave can be decomposed
http://arxiv.org/abs/0704.1407v1
into a linear combination of plane waves - it is a coher-
ent superposition of (an in principle infinite number of)
plane waves. The wave function of the fast electron can
be thus written as
ψ(r) =
ǫ(j)C(j)g e
i(k(j)+g)·r (3)
for incident wave and
ψ′(r) =
ǫ(l)D
(l)+h)·r (4)
for outgoing wave, where C
g , D
are so called Bloch co-
efficients, ǫ(j) (ǫ(l)) determine the excitation of the Bloch
wave with index j (l) and wave vector k(j) (k(l)) and g
(h) is a vector of the reciprocal lattice.
When we derive the Born approximation of DDSCS
starting with such fast electron wave functions, we will
obtain a sum of two kinds of terms: direct terms (DFFs)
as in the plane wave Born approximation Eq. (1), and
interference terms. These interference terms are a gener-
alization of the DFF - the mixed dynamic form factors11
(MDFFs). Each of them is defined by two wave vector
transfers, thus we label them S(q,q′, E).
The MDFF can be evaluated within a single particle
approximation as
S(q,q′, E) =
〈i|eiq·R̂|f〉〈f |e−iq
′·R̂|i〉δ(Ef − Ei − E)
where |i〉, |f〉 are the initial and final single-electron wave
functions of the target electron in the crystal. Thus
the definition of MDFF encompasses the notion of DFF,
Eq. (2), for q = q′. For more details about calculation
of MDFF see subsection II B.
The wave vector transfers are q
= k(l)−k(j)+h−g
and the total DDSCS will be a sum over all diads of q
and q′ vectors of terms
jlj′l′
ghg′h′
Sa(q,q
′, E)
q2q′2
where X
jlj′l′
ghg′h′
(a) is a the product of the coefficients of
the individual plane wave components of the fast electron
wave functions and a labels the position of the atoms
where the inelastic event can occur. The X
jlj′l′
ghg′h′
(a) co-
efficients are given by dynamical diffraction theory. This
will be covered in the next subsection IIA. The χf and
χ0 are magnitudes of wave vectors outside the crystal (in
the vacuum).
The calculation is thus split into two separate tasks.
i) Calculation of Bloch wave coefficients using the dy-
namical diffraction theory and identification of impor-
tant terms. This task is mainly ‘geometry dependent’,
although it can also contain some input from electronic
structure codes, namely the Coulomb part of crystal po-
tential. ii) Calculation of MDFFs requested by the dy-
namical diffraction theory. This part strongly depends
on the electronic structure of the studied system. The
final step is the summation of all terms.
A. Dynamical diffraction theory
The formalism, which will be described here is a gener-
alization of the formalism presented in Ref. 8,12 extend-
ing it beyond systematic row approximation by including
also higher-order Laue zones (HOLZ). The extension to
HOLZ is performed along lines presented in Refs. 13,14.
We will assume the high-energy Laue case, i.e. we can
safely neglect back-reflection and back-diffraction.
The Bloch wave vectors of the electron after entering
the crystal fulfill the continuity condition
(j) = χ+ γ(j)n (7)
where n is the unit vector normal to the crystal surface
and χ is the wave vector of the incoming electron. Only
the wave vector component normal to the surface can
change.
Expanding the wave function of the fast electron into
a linear combination of plane waves and substituting it
into the Schrödinger equation we obtain the secular equa-
tion14
K2 − (k(j) + g)2
h 6=0
 ei(k
(j)+g)·r = 0
where K2 = U0 + 2meE/~
2, m and e are, respectively,
the electron mass and charge, Ug = 2meVg/~
2 where
Vg are the Fourier components of the crystal potential,
which can be either calculated ab initio15,16 or obtained
from the tabulated forms of the potential17,18. It can
be shown14,19 that in the high energy limit the secular
equation, which is a quadratic eigenvalue problem in γ(j),
can be reduced to a linear eigenvalue problem AC(j) =
γ(j)C(j) where A is a non-hermitean matrix13,19
Agh =
K2 − (χ+ g)2
2(χ+ g) · n
δgh + (1− δgh)
2(χ+ g) · n
This eigenvalue problem can be transformed into a her-
mitean one using a diagonal matrix D with elements
Dgh = δgh
g · n
χ · n
Then the eigenvalue problem is equivalent to
(D1/2AD−1/2)(D1/2C(j)) = γ(j)(D1/2C(j)) or
(j) = γ(j)C̃(j), where the matrix Ã is hermitean
Ãgh =
K2 − (χ+ g)2
2(χ+ g) · n
δgh + (11)
+ (1− δgh)
[(χ+ g) · n][(χ+ h) · n]
and the original Bloch wave coefficients can be retrieved
using the relation
C(j)g = C̃
g · n
χ · n
By solving this eigenvalue problem we obtain the fast
electron wave function as a linear combination of eigen-
functions as given in Eq. (3). To obtain values for ǫ(j)
we need to impose boundary conditions, namely that the
electron is described by a single plane wave at the crys-
tal surface. The crystal surface is a plane defined by
the scalar product n · r = t0. Then the boundary con-
dition (in the high energy limit) leads to the following
condition14
ǫ(j) = C
−iγ(j)t0 (13)
It is easy to verify that
ψ(r)|n·r=t0 =
i(k(j)+g)·re−iγ
(j)t0
ei(χ+g)·r
n·re−iγ
(j)t0C
ei(χ+g)·r
ei(χ+g)·rδ0g
g · n
χ · n
= eiχ·r|n·r=t0 (14)
as required by the boundary condition. We have used
the continuity condition, Eq. (7), and the completeness
relation for the Bloch coefficients
δgh =
C̃(j)∗g C̃
= (15)
g · n
χ · n
h · n
χ · n
C(j)∗g C
Therefore the wave function of the fast moving electron
in the crystal, which becomes a single plane wave at n·r =
t0 is given by the following expression
ψ(r) =
iγ(j)(n·r−t0)ei(χ+g)·r (16)
The following discussion will be restricted to a partic-
ular case - a crystal with parallel surfaces. For such a
crystal with normals in the direction of the z axis we set
t0 = 0 for the fast electron entering the crystal and t0 = t
when leaving the crystal (t is the crystal thickness).
The inelastic event leads to a change of the energy
and momentum of the scattered electron. The detector
position determines the observed projection of the elec-
tron wave function (Bloch field) onto a plane wave af-
ter the inelastic event. Therefore the calculation of the
ELNES requires the solution of two independent eigen-
value problems describing an electron wave function be-
fore and after the inelastic event8,12. Invoking reciprocity
for electron propagation the outgoing wave can also be
considered as a time reversed solution of the Schrödinger
equation, also known as the reciprocal wave20 with the
source replacing the detector position.
Now we can identify the prefactors X
jlj′l′
ghg′h′
(a) from
Eq. (6). For the sake of clarity we will keep C
g for the
Bloch coefficients of the incoming electron and we use
for the Bloch coefficients of the outgoing electron en-
tering the detector (obtained from the two independent
eigenvalue problems). Similarly, superscript indices (j)
and (l) indicate eigenvalues and Bloch-vectors for incom-
ing and outgoing electron, respectively. We thus obtain
jlj′l′
ghg′h′
(a) = C
(j′)⋆
(l′)⋆
× ei(γ
(l)−γ(l
′))tei(q−q
′)·a (17)
where
q = k(l) − k(j) + h− g
′ = k(l
′) − k(j
′) + h′ − g′ (18)
In crystals the position of each atom can be decom-
posed into a sum of a lattice vector and a base vector,
a = R + u. Clearly, MDFF does not depend on R, but
only on u. It is then possible to perform analytically
the sum over all lattice vectors R under the approxima-
tion that the MDFF does not depend strongly on the j, l
indices. This is indeed a very good approximation, as
verified by numerical simulations (see below).
First we will treat the summation over all lattice vec-
tors. The sum in Eq. 6 can be separated into two terms
ei(q−q
′)·a =
ei(q−q
′)·u 1
ei(q−q
′)·R(19)
Since
q− q′ = [(γ(j) − γ(j
′))− (γ(l) − γ(l
′))]n
+ h− h′ + g′ − g (20)
and the algebric sum of g,h is simply a reciprocal lat-
tice vectors G, which fulfills eiG·R = 1, it is possible to
simplify the second term
ei(q−q
′)·R =
ei[(γ
(j)−γ(j
′))−(γ(l)−γ(l
′))]n·R (21)
For general orientations of the vector n this sum is
difficult to evaluate. In particular coordinate system with
n ‖ z and crystal axes a, b ⊥ z this sum leads8,21 to
ei(q−q
′)·R = NRe
i∆t/2 sin∆
so that the total sum over all atomic positions is
ei(q−q
′)·a = ei∆
sin∆ t
ei(q−q
′)·u(23)
where ∆ = (γ(j)−γ(j
′))− (γ(l)−γ(l
′)). The final expres-
sion of the DDSCS we write as
ghg′h′
Su(q,q
′, E)
q2q′2
ei(q−q
jlj′ l′
jlj′ l′
ghg′h′
Tjlj′l′(t) (24)
where
jlj′l′
ghg′h′
(j′)⋆
(l′)⋆
depends only on the eigenvectors of the incoming and
outgoing beam and
Tjlj′l′(t) = e
i[(γ(j)−γ(j
′))+(γ(l)−γ(l
′))] t
sin∆ t
is a thickness and eigenvalue dependent function.
Perturbative treatment of the absorption can be easily
introduced. If we denote by U ′g the absorptive part of the
potential, within the first order perturbation theory the
Bloch coefficients will not change, just the eigenvalues
will be shifted by iη(j) or iη(l) for the incoming or outgo-
ing wave, respectively. Particular η(j) can be calculated
using the following expression14
η(j) =
g,h U
g (χ+ g) · n
and similarly for the outgoing beam.
This way the eigenvalues change from γ(j) to γ(j)+iη(j)
and the ∆ acquires an imaginary part. Such approx-
imative treatment of absorption thus affects only the
thickness-dependent function Tjlj′l′(t).
Here we add a few practical considerations, which we
applied in our computer code. The sum in Eq. (24) is
performed over 8 indices for every energy and thickness
value. Such summation can easily grow to a huge number
of terms and go beyond the computational capability of
modern desktop computers. For example, if we assume
the splitting of the incoming (and outgoing) beam into
only 10 plane wave components, taking into account the
10 most strongly excited Bloch waves, we would have 108
terms per each energy and thickness. A calculation with
an energy mesh of 100 points at 100 different thicknesses
would include one trillion terms and require a consider-
able amount of computing time. However most of these
terms give a negligible contribution to the final sum.
Therefore several carefully chosen cut-off conditions are
required to keep the computing time reasonable without
any significant degradation of the accuracy.
The first cut-off condition used is based on the Ewald’s
sphere construction. Only plane wave components with
k + g close to the Ewald’s sphere will be excited. The
strength of the excitation decreases also with decreasing
crystal potential component Ug. A dimensionless param-
eter wg = sgξg - product of the excitation error and the
extinction distance14 - reflects both these criteria. There-
fore we can filter the list of beams by selecting only beams
with wg < wmax. Experience shows that in the final sum-
mation a fairly low number of beams is necessary to have
a well converged results (in systematic row conditions
this number is typically around 10). The convergence of
the corresponding Bloch coefficients requires solving an
eigenvalue problem with a much larger set of beams (sev-
eral hundreds). Therefore we defined two cut-off param-
eters for wg - the first for the solution of the eigenvalue
problem (typically wmax,1 is between 1000 and 5000) and
the second for the summation (wmax,2 typically between
50 and 100).
The second type of cut-off conditions is applied to
selection of Bloch waves, which enter the summation.
Once the set of beams for summation is determined, this
amounts to sorting the Bloch waves according to a prod-
uct of their excitation ǫ(j) and their norm on the subspace
defined by selected subset of beams, C
0 ||C
(j)||subsp. In
the systematic row conditions this value is large only for
a small number of Bloch waves. Typically in the exper-
imental geometries used for detection of EMCD one can
perform a summation over less than 10 Bloch waves to
have a well converged result (often 5 or 6 Bloch waves
are enough).
B. Mixed dynamic form factor
It can be seen from Eq. (5) that the calculation of the
MDFF requires the evaluation of two matrix elements be-
tween initial and final states of the target electron. The
derivation of the expression for the MDFF describing a
transition from core state nlκ (n, l, κ are the main, or-
bital and relativistic quantum numbers, respectively) to
a band state with energy E is presented in detail in the
supplementary material of Ref. 7 and in Ref. 8. Though,
note that in Ref. 8 the initial states are treated classi-
cally, which leads to somewhat different expression for
MDFF giving incorrect L2 − L3 branching ratio.
0 0.2 0.4 0.6 0.8 1
 ( a.u.
-0.04
-0.02
Re[ MDFF ] @ L
Re[ MDFF ] @ L
Im[ MDFF ] @ L
Im[ MDFF ] @ L
DFF @ L
DFF @ L
FIG. 1: Dependence of S(q, E) (top) and S(q,q′, E) with
q′ = G+q (bottom) on qz, calculated for the L2,3 edge of hcp-
Co, with G = (100), qx = −q
x = −|G|/2, qy = q
y = |G|/2.
The ratio between values calculated at L3 or L2 is constant
and equal to 2.1 for the real part and to −1 for the imaginary
part.
The final expression is7
S(q,q′, E) =
L′M ′S′
4πiλ−λ
(2l + 1)
[λ, λ′, L, L′]
× Y λµ (q/q)
µ′ (q
′/q′)〈jλ(q)〉ELSj〈jλ′ (q
′)〉EL′S′j
l λ L
0 0 0
l λ′ L′
0 0 0
l λ L
−m µ M
l λ′ L′
−m′ µ′ M ′
(−1)m+m
(2j + 1)
m S −jz
m′ S′ −jz
DLMS(νk)DL′M ′S′(νk)
∗δ(E + Enlκ − Eνk)
Here we made use of Wigner 3j-symbols, Y λµ are spheri-
cal harmonics, 〈jλ(q)〉ELSj are radial integrals of all the
radial-dependent terms (radial part of the wave function
of the core and band states, radial terms of the Rayleigh
expansion) and DLMS(νk) is the projection of the (νk)
Bloch state onto the LMS subspace within the atomic
sphere of the excited atom. For more details we refer to
the supplementary material of Ref. 7.
For evaluation of the radial integrals and Bloch state
projections DLMS(νk) we employ the density functional
theory22 within the local spin density approximation23.
-10 -5 0 5 10
 ( a.u.
-10 -5 0 5 10
 ( a.u.
Total
Dipole
λ = 0
λ = 1
λ = 2
a) b)
e) f)
FIG. 2: (online color) Decomposition of MDFF and dipole ap-
proximation calculated for hcp-Co with q′ − q = G = (100)
and qy = q
y = |G|/2 as a function of qx @ L3. Left column
- graphs a), c) and e), show S(q,q′, E) and right column,
graphs b), d) and f), show S(q,q′, E)/q2q′2. Top row - a)
and b) - is the DFF, middle row - c) and d) - is the real part
of MDFF and bottom row - e) and f) - is the imaginary part
of MDFF. The y-axes are in arbitrary units, but consistent
within the given column. The values for the L2 edge differ
only by a factor of 2.1 for the real part and −1 for the imag-
inary part. Note that the contributions of λ = 0 and 2 are
always negligible. See text for more details.
In Section IIA we used an approximation of negligi-
ble dependence of MDFF on the j, l indices (see Eq.20).
Generally, as the wave vector k(j,l) for each Bloch wave
changes slightly by an amount given by the corresponding
eigenvalue γ(j,l), the values of qz and q
z would change ac-
cordingly and therefore we should not be allowed to take
MDFF out of the sum over the indices j, l in the Eq. (24).
However, the change in qz (and q
z) induced by the eigen-
values γ(j,l) is small and can be neglected with respect
to the qz = χ0E/2E0 given by the energy loss E
26. To
demonstrate this we plot the dependence of MDFF on
qz, q
z for qx and qy corresponding to the main DFF and
MDFF terms, see Fig. 1. If qz is given in a.u.
−1 (atomic
units, 1 a.u.= 0.529178Å), typical values for L2,3 edges
of Fe, Co and Ni are around tenth of a.u.−1, whereas
typical values of γ(j,l) for strongly excited Bloch waves
are one or two orders of magnitude smaller. Thus the
approximation of weak j, l dependence of MDFF is well
justified.
Besides γ(j,l), the other factors determining the value
of qz are the energy of the edge, i.e. the energy lost by
the probe electron, the tilt with respect to the zone axis
and whether the excited beam is in a HOLZ. These last
factors have been included in our calculation. Only the
variations due to γ(j,l) are neglected, thus giving rise to
an error . 1%. If a more accurate treatment would be
needed, the smooth behavior of MDFF with respect to
qz would allow to use simple linear or quadratic interpo-
lation/extrapolation methods.
As mentioned in the introduction and explained in
Refs. 1,7, dichroism in the TEM is made possible by
the analogous role that the polarization vector ε and the
wave vector transfer q play in the dipole approximation
of the DDSCS. However we do not restrict our calcu-
lations to the dipole approximation. We use the more
complete expression Eq. (5).
To evaluate the accuracy of the dipole approximation,
we compare the dipole approximation of MDFF with the
full calculation (with λ up to 3) also showing λ-diagonal
components of the MDFF, Fig. 2. Because the domi-
nant contribution to the signal originates from (dipole
allowed) 2p→ 3d transitions, the λ = λ′ = 1 term nearly
coincides with the total MDFF. While the dipole approx-
imation works relatively well for the studied systems,
particularly the MDFF divided by squares of momen-
tum transfer vectors (right column of the Fig. 2), it has
significantly different asymptotic behaviours for larger q-
vectors. The λ = λ′ = 1 term provides a much better
approximation, which remains very accurate also in the
large q region.
It is worth mentioning that thanks to the properties
of the Gaunt coefficients the 2p → 3d transitions are all
included in the λ = 1 and λ = 3 contributions. Thanks
to the negligible value of the radial integrals for λ =
3 the terms with λ = 1 account for the large majority
of the calculated signal. The contributions from λ =
0, 2 describe transitions from 2p to valence p or f states
and are always negligible due to the composition of the
density of states beyond the Fermi level. They practically
overlap with the zero axis in all the six parts of Fig. 2.
It can be shown24 that in the dipole approximation
the real part of the MDFF is proportional to q · q′ and
the imaginary part is proportional to q × q′. A little
algebra can thus show that the imaginary part of the
MDFF is, in the geometry described in the caption of
Fig. 2, constant with respect to qx. As expected, the DFF
(which is proportional to q2) has a minimum at qx = 0,
where S(q, E)/q4 has a maximum. For the MDFF (and
corresponding S(q,q′, E)/q2q′2) these minima and max-
ima are centered at qx = −G/2 = −0.76 a.u.
−1 where
|qx| = |q
III. RESULTS
We summarize the results obtained for body-centered
cubic iron (bcc-Fe), hexagonal close-packed cobalt (hcp-
Co) and face-centered cubic nickel (fcc-Ni) crystals,
which are also the first samples prepared for EMCD mea-
surements. These results are valuable for optimization of
the experimental setup.
The geometry setup for observing the dichroic effect7
consists in creating a two-beam case by tilting the beam
away from a zone axis (here (001)) by a few degrees and
then setting the Laue circle center equal to G/2 for the
G vector to be excited. In analogy to XMCD, where two
measurements are performed for left- and righ-handed
circularly polarized light, here we perform two measure-
ments by changing the position of the detector, which
lies once at the top and once at the bottom of the Thales
circle having as diameter the line connecting the diffrac-
tion spots 0 and G. This geometry setup, together with
the crystal structure, is an input for the calculation of
the Bloch wave coefficients (within the systematic row
approximation) using the dynamical diffraction theory
code described in section IIA.
The electronic structure was calculated using the
WIEN2k package15, which is a state-of-the-art implemen-
tation of the full-potential linearized augmented plane
waves method. The experimental values of lattice param-
eters were used. More than 10000 k-points were used to
achieve a very good converge of the Brillouin zone inte-
grations. Atomic sphere sizes were 2.2, 2.3 and 2.2 bohr
radii for bcc-Fe, hcp-Co and fcc-Ni, respectively. The re-
sulting electronic structure was the input for the calcula-
tion of the individual MDFFs required for the summation
(see Section II B).
In the three studied cases the dichroic effect is domi-
nated by the transitions to the unoccupied 3d states. The
d-resolved spin-up density of states (DOS) is almost fully
occupied, while the spin-down d-DOS is partially unoccu-
pied. In Fig. 3 we compare the d-DOS with the dichroic
signal at the L3 edge. Due to negligible orbital moments
in these compounds the L2 edge shows a dichroic signal
of practically the same magnitude but with opposite sign.
The shape of the calculated dichroic peaks corresponds to
the difference of spin-up and spin-down d-DOS, similarly
to XMCD, as it was shown for the same set of systems in
Ref. 25. The calculations were performed within system-
atic row conditions with G = (200) for bcc-Fe and fcc-Ni
and G = (100) for hcp-Co. The sample thicknesses were
set to 20 nm, 10 nm and 8 nm for bcc-Fe, hcp-Co and fcc-
Ni, respectively. These values were found to be optimal
for these systems in the given experimental geometry.
An interesting point is the comparison of the strength
of the dichroic signal. According to the d-DOS projec-
tions one would expect comparable strength of signals
for the three elements under study. But the dichroic
signal of hcp-Co seems to be approximately a factor of
two smaller than that of the other two. The reason for
that can be explained by simple geometrical considera-
-6 -4 -2 0 2 4
E - E
 ( eV )
0 2 4
 ( eV )
FIG. 3: (online color) Spin-resolved d-densities of states (left) and resulting signal on L3 edge (right) on bcc-Fe, hcp-Co and
fcc-Ni (from top to bottom) at optimal thickness (see text). Spin-up DOS is drawn using a solid black line (positive) and spin-
down DOS using a dashed red line (negative). DDSCS for the (+) detector position is drawn using a solid blue line, DDSCS
for the (-) position is drawn using a dashed green line. The dichroic signal (difference) is the hatched red area. G = (200) for
bcc-Fe and fcc-Ni and (100) for hcp-Co.
tions starting from Eq. (24). For simplicity we consider
only the main contributions: the DFF S(q,q, E) and
the MDFF S(q,q′, E) with q ⊥ q′. For bcc-Fe and fcc-
Ni the summation over u within the Bravais cell leads
always to the structure factor 2 and 4, respectively, be-
cause q′ − q = G is a kinematically allowed reflection.
This factor cancels out after division by the number of
atoms in the Bravais cell. Therefore it does not matter,
what is the value of q-vectors, the sum over the atoms
is equal to S(q,q′, E)/q2q′2 itself. On the other hand,
the unit cell of hcp-Co contains two equivalent atoms at
positions u1 = (
) and u2 = (
). For the two
DFFs q = q′ and the exponential reduces to 1; since there
are two such terms, after division by Nu the sum equals
again the DFF itself. But for the main MDFF we have
q ⊥ q′ and the exponential factor will in general weight
the terms. One can easily see, that q′ − q = G. For the
G = (100) systematic row case, which was used for cal-
culation of hcp-Co in Fig. 3 the exponentials evaluate to
the complex numbers − 1
and − 1
for u1 and
u2, respectively. Because of symmetry, the MDFFs for
both atoms are equal and then the sum 1
leads to a
factor − 1
for the MDFF contribution, i.e. the influence
of its imaginary part, which is responsible for dichroism,
on the DDSCS is reduced by a factor of two.
To optimize the dichroic signal strength of hcp-Co, we
require G · u1 = G · u2 = 2πn, which gives in principle
an infinite set of possible G vectors. The one with lowest
hkl indices is G = (110). A calculation for this geometry
setup leads to approximately twice the dichroic signal,
see Fig. 4 and compare to the corresponding graph in
Fig. 3.
For the optimization of the experimental setup it is
important to know how sensitive the results are to varia-
tion of the parameters like the thickness of the sample or
the accuracy of the detector position. Another question
related to this is also the sensitivity to the finite size of
the convergence and collection angles α and β. In the
following text we will address these questions.
The thickness influences the factor Tjlj′l′ in the
Eq. (24) only. This factor leads to the so called pen-
dellösung oscillations - modulations of the signal strength
-1 0 1 2 3 4
Energy loss ( eV )
dichroism (%)
DDSCS(+)
DDSCS(-)
FIG. 4: L3 peak of hcp-Co calculated for the G = (110)
systematic row at 18nm. See caption of Fig. 3. The peaks
have been renormalized so that their sum is 1, therefore their
difference is the dichroic signal (ca. 15% in this case).
as a function of thickness. This also influences the
strength of the dichroic signal. Results of such calcula-
tions are displayed in Fig. 5 (we did not include absorp-
tion into these simulations, so that all signal variations
are only due to the geometry of the sample). From these
simulations it follows that a well defined thickness of the
sample is a very important factor. Relatively small vari-
ations of the thickness can induce large changes in the
dichroic signal, particularly in fcc-Ni. From the figure
one can deduce that the optimal thickness for a bcc-Fe
sample should be between 8 nm and 22 nm (of course, due
to absorption, thinner samples within this range would
have a stronger signal), for hcp-Co between 15 nm and
22 nm and for fcc-Ni it is a relatively narrow interval -
between 6 nm and 10 nm. However, we stress that these
results depend on the choice of the systematic row vector
G. For example hcp-Co with G = (100) (instead of (110)
shown in Fig. 5) has a maximum between 5 nm and 15
nm (although it is much lower, as discussed before).
Taking the optimal thickness, namely 20 nm, 18 nm
and 8 nm for bcc-Fe, hcp-Co and fcc-Ni, respectively, we
calculated the dependence of the dichroic signal on the
detector position. We particularly tested changes of the
dichroic signal when the detector is moved away from its
default position in the direction perpendicular to G, see
Fig. 6. It is interesting to note that the maximum abso-
lute difference occurs for a value of qy smaller than |G|/2.
This can be qualitatively explained by considering the
non-zero value of qz and q
z, i.e. q and q
′ are not exactly
perpendicular at the default detector positions. More-
over the MDFF enters the summation always divided by
q2q′2 and the lengths of q-vectors decrease with decreas-
ing qy. The important message we can deduce from this
figure is that the dichroic signal is only weakly sensitive
to the accuracy of qy since even displacement by 10-20%
in the detector default qy positions (i.e. qy = ±G/2) do
not affect significantly the measured dichroic signal.
0 20 40 60 80 100
Thickness ( nm )
FIG. 5: (online color) Dependence of the DDSCS and of the
dichroic signal on sample thickness for a) bcc-Fe, b) hcp-Co
and c) fcc-Ni. Systematic row vector G = (200) was used for
bcc-Fe and fcc-Ni, while for hcp-Co G = (110) was chosen.
The blue and green solid curves are DDSCSs calculated for
the (+) and (-) detector positions, the dashed black curve is
the DFF part of the DDSCS (it is identical for both detector
positions). The red line with circles is the relative dichroism
defined as difference of DDSCSs divided by their sum, the red
solid curve is the absolute dichroism - difference of DDSCSs.
Related to this is a study of the dependence of the
dichroic signal on the finite size of the convergence and
collection angles α and β. We performed a calculation
for the three studied metals and found that collection
and convergence half-angles up to 2 mrad weakens the
relative dichroic signal by less than 10%.
IV. CONCLUSIONS
We have developed a computer code package for the
calculation of electron energy loss near edge spectra,
which includes the theory of dynamical Bragg diffrac-
tion. We applied the code to the recently discovered phe-
nomenon of magnetic chiral dichroism in the TEM and
we demonstrated the relation of the dichroic peak shape
to the difference of d-projections of the spin-resolved
density of states in analogy with similar observation for
XMCD.
Using this code we examined the validity of the dipole
approximation, which is often assumed. We found that
-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8
 / |G| ( dimensionless )
bcc-Fe
hcp-Co
fcc-Ni
absolute dichroism (arb. units)
FIG. 6: (online color) Dependence of the dichroic signal on
detector displacements along qy . The full symbols correspond
to the relative dichroic signal while the open symbols to the
difference of the DDSCS for both detector positions. These
are in arbitrary units and their magnitudes are not directly
comparable. Vertical lines are showing the default detector
positions.
for the 3d ferromagnetic systems studied it is a reasonable
approximation, however with wrong asymptotic proper-
ties - it overestimates the contributions from larger q-
vectors. A very accurate approximation for the studied
systems is the λ = λ′ = 1 approximation, which treats
appropriately the dominant p→ d dipole transitions and
remains very accurate also for large q, q′.
In order to provide guidance to the experimentalist we
have investigated the strength of the dichroic signal as a
function of the sample thickness and the precision of the
detector placement. While the dichroic signal strength is
rather robust with respect to the precision of the detector
placement, the thickness of the specimen influences the
signal considerably. Therefore it might be a challenge to
produce samples with optimum thickness and selecting
the best systematic row Bragg spot. Our calculations
yield best thicknesses in order to detect EMCD of the
iron and nickel samples for the systematic rowG = (200)
to be 8-22 nm and 6-10 nm, respectively, and for cobalt
in the systematic row G = (110) to be 15-22 nm.
Acknowledgments
We thank Dr. Cécile Hébert and Dr. Pavel Novák for
stimulating discussions. This work has been supported
by the European Commission, contract nr. 508971 (CHI-
RALTEM).
∗ Electronic address: rusz@mag.mff.cuni.cz
† currently at Department of Physics, Uppsala University,
Box 530, S-751 21 Uppsala, Sweden
1 A. P. Hitchcock, Jpn. J. Appl. Phys. 32(2), 176 (1993).
2 J. Yuan and N. K. Menon, J. Appl. Phys. 81(8), 5087
(1997).
3 B. T. Thole, P. Carra, F. Sette and G. van der Laan, Phys.
Rev. Lett. 68, 1943 (1992).
4 S. W. Lovesey and S. P. Collins, X-Ray Scattering and Ab-
sorption by Magnetic Materials, Clarendon Press, Oxford,
UK, 1996.
5 J. Stöhr, J. Electron Spectrosc. Relat. Phenom. 75, 253
(1995).
6 C. Hébert and P. Schattschneider, Ultramicroscopy 96, 463
(2003).
7 P. Schattschneider, S. Rubino, C. Hébert, J. Rusz, J.
Kuneš, P. Novák, E. Carlino, M. Fabrizioli, G. Panaccione
and G. Rossi, Nature 441, 486 (2006).
8 M. Nelhiebel, Ph.D. Thesis, 1999, Vienna University of
Technology.
9 M. Nelhiebel, P. Schattschneider and B. Jouffrey, Phys.
Rev. Lett. 85(9), 1847 (2000).
10 P. Schattschneider, C. Hébert, H. Franco and B. Jouffrey,
Phys. Rev. B 72, 045142 (2005).
11 H. Kohl and H. Rose, Advances in electronics and electron
optics 65, 173 (1985).
12 P. Schattschneider, M. Nelhiebel and B. Jouffrey, Phys.
Rev. B 59, 10959 (1999).
13 L.-M. Peng and M. J. Whelan, Proc. R. Soc. London Ser.
A 431, 111 (1990).
14 A. J. F. Metherell, in Electron Microscopy in Materials
Science II, 397 (1975), CEC, Luxembourg.
15 P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka and
J. Luitz, 2001 WIEN2k, Vienna University of Technology
(ISBN 3-9501031-1-2).
16 The WIEN2k package can calculate X-ray structure fac-
tors. By supplying the potential instead of charge density
it is possible to use the same code to calculate the electron
structure factors using a simple script.
17 A. Weickenmeier and H. Kohl, Acta Cryst. A47, 590
(1991).
18 P. A. Doyle and P. S. Turner, Acta Cryst. A24, 390 (1968).
19 A. L. Lewis, R. E. Villagrana and A. J. F. Metherell, Acta
Cryst. A34, 138 (1978).
20 Y. Kainuma, Acta Cryst. 8, 247 (1955).
21 P. Schattschneider and W. S. M. Werner, J. Electron Spec-
trosc. Relat. Phenom. 143, 81 (2005).
22 P. Hohenberg, W. Kohn, Phys. Rev. 136, B864 (1964) and
W. Kohn, L. J. Sham, Phys. Rev. 140, A1133 (1965).
23 J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
24 P. Schattschneider, C. Hébert, S. Rubino, M. Stöger-
Pollach, J. Rusz, P. Novák, submitted.
25 J. Kuneš and P. M. Oppeneer, Phys. Rev. B 67, 024431
(2003).
26 The question of momentum conservation in the z direction
in the inelastic interaction in a crystal of finite thickness is
related to the probability of inter- and intrabranch transi-
tions of the probe electron?
mailto:rusz@mag.mff.cuni.cz

  Recently it was demonstrated (Schattschneider et al., Nature 441 (2006),
486), that an analogue of the X-ray magnetic circular dichroism (XMCD)
experiment can be performed with the transmission electron microscope (TEM).
The new phenomenon has been named energy-loss magnetic chiral dichroism (EMCD).
In this work we present a detailed ab initio study of the chiral dichroism in
the Fe, Co and Ni transition elements. We discuss the methods used for the
simulations together with the validity and accuracy of the treatment, which
can, in principle, apply to any given crystalline specimen. The dependence of
the dichroic signal on the sample thickness, accuracy of the detector position
and the size of convergence and collection angles is calculated.

First principles theory of chiral dichroism in electron microscopy applied to 3d
ferromagnets
Ján Rusz∗
Institute of Physics, Academy of Sciences of the Czech Republic,
Cukrovarnická 10, 162 53 Prague 6, Czech Republic†
Stefano Rubino and Peter Schattschneider
Institute for Solid State Physics, Vienna University of Technology,
Wiedner Hauptstrasse 8-10/138, A-1040 Vienna, Austria
(Dated: August 2, 2021)
Recently it was demonstrated (Schattschneider et al., Nature 441 (2006), 486), that an analogue of
the X-ray magnetic circular dichroism (XMCD) experiment can be performed with the transmission
electron microscope (TEM). The new phenomenon has been named energy-loss magnetic chiral
dichroism (EMCD). In this work we present a detailed ab initio study of the chiral dichroism in the
Fe, Co and Ni transition elements. We discuss the methods used for the simulations together with
the validity and accuracy of the treatment, which can, in principle, apply to any given crystalline
specimen. The dependence of the dichroic signal on the sample thickness, accuracy of the detector
position and the size of convergence and collection angles is calculated.
PACS numbers:
Keywords: density functional theory, chiral dichroism, transmission electron microscopy, dynamical diffrac-
tion theory
I. INTRODUCTION
The analogy between X-ray absorption spectroscopy
(XAS) and electron energy loss spectroscopy (EELS) has
been recognized long ago1,2. The role of the polarization
vector ε in XAS is similar to the role of the wave vector
transfer q in EELS. This has made feasible the detection
of linear dichroism in the TEM. However the counter-
part of X-ray magnetic circular dichroism (XMCD)3,4,5
experiments with electron probes was thought to be tech-
nically impossible due to the low intensity of existing
spin polarized electron sources. XMCD is an important
technique providing atom-specific information about the
magnetic properties of materials. Particularly the near
edge spectra, where a well localized strongly bound elec-
tron with l 6= 0 is excited to an unoccupied band state,
allow to measure spin and orbital moments. Soon after
the proposal of an experimental setup for detection of cir-
cular dichroism using a standard non-polarized electron
beam in the TEM6 it was demonstrated that such ex-
periments (called energy-loss magnetic chiral dichroism,
EMCD) are indeed possible7. This novel technique is of
considerable interest for nanomagnetism and spintron-
ics according to the high spatial resolution of the TEM.
However, its optimization involves many open questions.
In this work we provide theoretical ab initio predictions
of the dependence of the dichroic signal in the EMCD
experiment on several experimental conditions, such as
sample thickness, detector placement and the finite size
of convergence and collection angles. This information
should help to optimize the experimental geometry in
order to maximize the signal to noise ratio.
The structure of this work is as follows: In section II
we first describe the computational approach based on
the dynamical diffraction theory and electronic structure
calculations. We also discuss the validity of several ap-
proximations for the mixed dynamic form factor. In sec-
tion III we study the dependence of the dichroic signal of
bcc-Fe, hcp-Co and fcc-Ni on various experimental con-
ditions. This section is followed by a concluding section
summarizing the most important findings.
II. METHOD OF CALCULATION
We will follow the derivations of the double differential
scattering cross-section (DDSCS) presented in Refs. 8,9.
Within the first-order Born approximation10 DDSCS is
written as
S(q, E)
S(q, E) =
|〈i|eiq·R̂|f〉|2δ(Ef − Ei − E) (2)
where q = kf−k0 is the difference (wave vector transfer)
between final wave vector kf and initial wave vector k0
of the fast electron; γ = 1/
1− v2/c2 is a relativistic
factor and a0 is Bohr radius. The S(q, E) is the so called
dynamic form factor (DFF)11.
This equation is valid only if the initial and final wave
functions of the fast electron are plane waves. In the
crystal the full translation symmetry is broken and as
a result, the electron wave function becomes a superpo-
sition of Bloch waves, which reflects the discrete trans-
lation symmetry. Each Bloch wave can be decomposed
http://arxiv.org/abs/0704.1407v1
into a linear combination of plane waves - it is a coher-
ent superposition of (an in principle infinite number of)
plane waves. The wave function of the fast electron can
be thus written as
ψ(r) =
ǫ(j)C(j)g e
i(k(j)+g)·r (3)
for incident wave and
ψ′(r) =
ǫ(l)D
(l)+h)·r (4)
for outgoing wave, where C
g , D
are so called Bloch co-
efficients, ǫ(j) (ǫ(l)) determine the excitation of the Bloch
wave with index j (l) and wave vector k(j) (k(l)) and g
(h) is a vector of the reciprocal lattice.
When we derive the Born approximation of DDSCS
starting with such fast electron wave functions, we will
obtain a sum of two kinds of terms: direct terms (DFFs)
as in the plane wave Born approximation Eq. (1), and
interference terms. These interference terms are a gener-
alization of the DFF - the mixed dynamic form factors11
(MDFFs). Each of them is defined by two wave vector
transfers, thus we label them S(q,q′, E).
The MDFF can be evaluated within a single particle
approximation as
S(q,q′, E) =
〈i|eiq·R̂|f〉〈f |e−iq
′·R̂|i〉δ(Ef − Ei − E)
where |i〉, |f〉 are the initial and final single-electron wave
functions of the target electron in the crystal. Thus
the definition of MDFF encompasses the notion of DFF,
Eq. (2), for q = q′. For more details about calculation
of MDFF see subsection II B.
The wave vector transfers are q
= k(l)−k(j)+h−g
and the total DDSCS will be a sum over all diads of q
and q′ vectors of terms
jlj′l′
ghg′h′
Sa(q,q
′, E)
q2q′2
where X
jlj′l′
ghg′h′
(a) is a the product of the coefficients of
the individual plane wave components of the fast electron
wave functions and a labels the position of the atoms
where the inelastic event can occur. The X
jlj′l′
ghg′h′
(a) co-
efficients are given by dynamical diffraction theory. This
will be covered in the next subsection IIA. The χf and
χ0 are magnitudes of wave vectors outside the crystal (in
the vacuum).
The calculation is thus split into two separate tasks.
i) Calculation of Bloch wave coefficients using the dy-
namical diffraction theory and identification of impor-
tant terms. This task is mainly ‘geometry dependent’,
although it can also contain some input from electronic
structure codes, namely the Coulomb part of crystal po-
tential. ii) Calculation of MDFFs requested by the dy-
namical diffraction theory. This part strongly depends
on the electronic structure of the studied system. The
final step is the summation of all terms.
A. Dynamical diffraction theory
The formalism, which will be described here is a gener-
alization of the formalism presented in Ref. 8,12 extend-
ing it beyond systematic row approximation by including
also higher-order Laue zones (HOLZ). The extension to
HOLZ is performed along lines presented in Refs. 13,14.
We will assume the high-energy Laue case, i.e. we can
safely neglect back-reflection and back-diffraction.
The Bloch wave vectors of the electron after entering
the crystal fulfill the continuity condition
(j) = χ+ γ(j)n (7)
where n is the unit vector normal to the crystal surface
and χ is the wave vector of the incoming electron. Only
the wave vector component normal to the surface can
change.
Expanding the wave function of the fast electron into
a linear combination of plane waves and substituting it
into the Schrödinger equation we obtain the secular equa-
tion14
K2 − (k(j) + g)2
h 6=0
 ei(k
(j)+g)·r = 0
where K2 = U0 + 2meE/~
2, m and e are, respectively,
the electron mass and charge, Ug = 2meVg/~
2 where
Vg are the Fourier components of the crystal potential,
which can be either calculated ab initio15,16 or obtained
from the tabulated forms of the potential17,18. It can
be shown14,19 that in the high energy limit the secular
equation, which is a quadratic eigenvalue problem in γ(j),
can be reduced to a linear eigenvalue problem AC(j) =
γ(j)C(j) where A is a non-hermitean matrix13,19
Agh =
K2 − (χ+ g)2
2(χ+ g) · n
δgh + (1− δgh)
2(χ+ g) · n
This eigenvalue problem can be transformed into a her-
mitean one using a diagonal matrix D with elements
Dgh = δgh
g · n
χ · n
Then the eigenvalue problem is equivalent to
(D1/2AD−1/2)(D1/2C(j)) = γ(j)(D1/2C(j)) or
(j) = γ(j)C̃(j), where the matrix Ã is hermitean
Ãgh =
K2 − (χ+ g)2
2(χ+ g) · n
δgh + (11)
+ (1− δgh)
[(χ+ g) · n][(χ+ h) · n]
and the original Bloch wave coefficients can be retrieved
using the relation
C(j)g = C̃
g · n
χ · n
By solving this eigenvalue problem we obtain the fast
electron wave function as a linear combination of eigen-
functions as given in Eq. (3). To obtain values for ǫ(j)
we need to impose boundary conditions, namely that the
electron is described by a single plane wave at the crys-
tal surface. The crystal surface is a plane defined by
the scalar product n · r = t0. Then the boundary con-
dition (in the high energy limit) leads to the following
condition14
ǫ(j) = C
−iγ(j)t0 (13)
It is easy to verify that
ψ(r)|n·r=t0 =
i(k(j)+g)·re−iγ
(j)t0
ei(χ+g)·r
n·re−iγ
(j)t0C
ei(χ+g)·r
ei(χ+g)·rδ0g
g · n
χ · n
= eiχ·r|n·r=t0 (14)
as required by the boundary condition. We have used
the continuity condition, Eq. (7), and the completeness
relation for the Bloch coefficients
δgh =
C̃(j)∗g C̃
= (15)
g · n
χ · n
h · n
χ · n
C(j)∗g C
Therefore the wave function of the fast moving electron
in the crystal, which becomes a single plane wave at n·r =
t0 is given by the following expression
ψ(r) =
iγ(j)(n·r−t0)ei(χ+g)·r (16)
The following discussion will be restricted to a partic-
ular case - a crystal with parallel surfaces. For such a
crystal with normals in the direction of the z axis we set
t0 = 0 for the fast electron entering the crystal and t0 = t
when leaving the crystal (t is the crystal thickness).
The inelastic event leads to a change of the energy
and momentum of the scattered electron. The detector
position determines the observed projection of the elec-
tron wave function (Bloch field) onto a plane wave af-
ter the inelastic event. Therefore the calculation of the
ELNES requires the solution of two independent eigen-
value problems describing an electron wave function be-
fore and after the inelastic event8,12. Invoking reciprocity
for electron propagation the outgoing wave can also be
considered as a time reversed solution of the Schrödinger
equation, also known as the reciprocal wave20 with the
source replacing the detector position.
Now we can identify the prefactors X
jlj′l′
ghg′h′
(a) from
Eq. (6). For the sake of clarity we will keep C
g for the
Bloch coefficients of the incoming electron and we use
for the Bloch coefficients of the outgoing electron en-
tering the detector (obtained from the two independent
eigenvalue problems). Similarly, superscript indices (j)
and (l) indicate eigenvalues and Bloch-vectors for incom-
ing and outgoing electron, respectively. We thus obtain
jlj′l′
ghg′h′
(a) = C
(j′)⋆
(l′)⋆
× ei(γ
(l)−γ(l
′))tei(q−q
′)·a (17)
where
q = k(l) − k(j) + h− g
′ = k(l
′) − k(j
′) + h′ − g′ (18)
In crystals the position of each atom can be decom-
posed into a sum of a lattice vector and a base vector,
a = R + u. Clearly, MDFF does not depend on R, but
only on u. It is then possible to perform analytically
the sum over all lattice vectors R under the approxima-
tion that the MDFF does not depend strongly on the j, l
indices. This is indeed a very good approximation, as
verified by numerical simulations (see below).
First we will treat the summation over all lattice vec-
tors. The sum in Eq. 6 can be separated into two terms
ei(q−q
′)·a =
ei(q−q
′)·u 1
ei(q−q
′)·R(19)
Since
q− q′ = [(γ(j) − γ(j
′))− (γ(l) − γ(l
′))]n
+ h− h′ + g′ − g (20)
and the algebric sum of g,h is simply a reciprocal lat-
tice vectors G, which fulfills eiG·R = 1, it is possible to
simplify the second term
ei(q−q
′)·R =
ei[(γ
(j)−γ(j
′))−(γ(l)−γ(l
′))]n·R (21)
For general orientations of the vector n this sum is
difficult to evaluate. In particular coordinate system with
n ‖ z and crystal axes a, b ⊥ z this sum leads8,21 to
ei(q−q
′)·R = NRe
i∆t/2 sin∆
so that the total sum over all atomic positions is
ei(q−q
′)·a = ei∆
sin∆ t
ei(q−q
′)·u(23)
where ∆ = (γ(j)−γ(j
′))− (γ(l)−γ(l
′)). The final expres-
sion of the DDSCS we write as
ghg′h′
Su(q,q
′, E)
q2q′2
ei(q−q
jlj′ l′
jlj′ l′
ghg′h′
Tjlj′l′(t) (24)
where
jlj′l′
ghg′h′
(j′)⋆
(l′)⋆
depends only on the eigenvectors of the incoming and
outgoing beam and
Tjlj′l′(t) = e
i[(γ(j)−γ(j
′))+(γ(l)−γ(l
′))] t
sin∆ t
is a thickness and eigenvalue dependent function.
Perturbative treatment of the absorption can be easily
introduced. If we denote by U ′g the absorptive part of the
potential, within the first order perturbation theory the
Bloch coefficients will not change, just the eigenvalues
will be shifted by iη(j) or iη(l) for the incoming or outgo-
ing wave, respectively. Particular η(j) can be calculated
using the following expression14
η(j) =
g,h U
g (χ+ g) · n
and similarly for the outgoing beam.
This way the eigenvalues change from γ(j) to γ(j)+iη(j)
and the ∆ acquires an imaginary part. Such approx-
imative treatment of absorption thus affects only the
thickness-dependent function Tjlj′l′(t).
Here we add a few practical considerations, which we
applied in our computer code. The sum in Eq. (24) is
performed over 8 indices for every energy and thickness
value. Such summation can easily grow to a huge number
of terms and go beyond the computational capability of
modern desktop computers. For example, if we assume
the splitting of the incoming (and outgoing) beam into
only 10 plane wave components, taking into account the
10 most strongly excited Bloch waves, we would have 108
terms per each energy and thickness. A calculation with
an energy mesh of 100 points at 100 different thicknesses
would include one trillion terms and require a consider-
able amount of computing time. However most of these
terms give a negligible contribution to the final sum.
Therefore several carefully chosen cut-off conditions are
required to keep the computing time reasonable without
any significant degradation of the accuracy.
The first cut-off condition used is based on the Ewald’s
sphere construction. Only plane wave components with
k + g close to the Ewald’s sphere will be excited. The
strength of the excitation decreases also with decreasing
crystal potential component Ug. A dimensionless param-
eter wg = sgξg - product of the excitation error and the
extinction distance14 - reflects both these criteria. There-
fore we can filter the list of beams by selecting only beams
with wg < wmax. Experience shows that in the final sum-
mation a fairly low number of beams is necessary to have
a well converged results (in systematic row conditions
this number is typically around 10). The convergence of
the corresponding Bloch coefficients requires solving an
eigenvalue problem with a much larger set of beams (sev-
eral hundreds). Therefore we defined two cut-off param-
eters for wg - the first for the solution of the eigenvalue
problem (typically wmax,1 is between 1000 and 5000) and
the second for the summation (wmax,2 typically between
50 and 100).
The second type of cut-off conditions is applied to
selection of Bloch waves, which enter the summation.
Once the set of beams for summation is determined, this
amounts to sorting the Bloch waves according to a prod-
uct of their excitation ǫ(j) and their norm on the subspace
defined by selected subset of beams, C
0 ||C
(j)||subsp. In
the systematic row conditions this value is large only for
a small number of Bloch waves. Typically in the exper-
imental geometries used for detection of EMCD one can
perform a summation over less than 10 Bloch waves to
have a well converged result (often 5 or 6 Bloch waves
are enough).
B. Mixed dynamic form factor
It can be seen from Eq. (5) that the calculation of the
MDFF requires the evaluation of two matrix elements be-
tween initial and final states of the target electron. The
derivation of the expression for the MDFF describing a
transition from core state nlκ (n, l, κ are the main, or-
bital and relativistic quantum numbers, respectively) to
a band state with energy E is presented in detail in the
supplementary material of Ref. 7 and in Ref. 8. Though,
note that in Ref. 8 the initial states are treated classi-
cally, which leads to somewhat different expression for
MDFF giving incorrect L2 − L3 branching ratio.
0 0.2 0.4 0.6 0.8 1
 ( a.u.
-0.04
-0.02
Re[ MDFF ] @ L
Re[ MDFF ] @ L
Im[ MDFF ] @ L
Im[ MDFF ] @ L
DFF @ L
DFF @ L
FIG. 1: Dependence of S(q, E) (top) and S(q,q′, E) with
q′ = G+q (bottom) on qz, calculated for the L2,3 edge of hcp-
Co, with G = (100), qx = −q
x = −|G|/2, qy = q
y = |G|/2.
The ratio between values calculated at L3 or L2 is constant
and equal to 2.1 for the real part and to −1 for the imaginary
part.
The final expression is7
S(q,q′, E) =
L′M ′S′
4πiλ−λ
(2l + 1)
[λ, λ′, L, L′]
× Y λµ (q/q)
µ′ (q
′/q′)〈jλ(q)〉ELSj〈jλ′ (q
′)〉EL′S′j
l λ L
0 0 0
l λ′ L′
0 0 0
l λ L
−m µ M
l λ′ L′
−m′ µ′ M ′
(−1)m+m
(2j + 1)
m S −jz
m′ S′ −jz
DLMS(νk)DL′M ′S′(νk)
∗δ(E + Enlκ − Eνk)
Here we made use of Wigner 3j-symbols, Y λµ are spheri-
cal harmonics, 〈jλ(q)〉ELSj are radial integrals of all the
radial-dependent terms (radial part of the wave function
of the core and band states, radial terms of the Rayleigh
expansion) and DLMS(νk) is the projection of the (νk)
Bloch state onto the LMS subspace within the atomic
sphere of the excited atom. For more details we refer to
the supplementary material of Ref. 7.
For evaluation of the radial integrals and Bloch state
projections DLMS(νk) we employ the density functional
theory22 within the local spin density approximation23.
-10 -5 0 5 10
 ( a.u.
-10 -5 0 5 10
 ( a.u.
Total
Dipole
λ = 0
λ = 1
λ = 2
a) b)
e) f)
FIG. 2: (online color) Decomposition of MDFF and dipole ap-
proximation calculated for hcp-Co with q′ − q = G = (100)
and qy = q
y = |G|/2 as a function of qx @ L3. Left column
- graphs a), c) and e), show S(q,q′, E) and right column,
graphs b), d) and f), show S(q,q′, E)/q2q′2. Top row - a)
and b) - is the DFF, middle row - c) and d) - is the real part
of MDFF and bottom row - e) and f) - is the imaginary part
of MDFF. The y-axes are in arbitrary units, but consistent
within the given column. The values for the L2 edge differ
only by a factor of 2.1 for the real part and −1 for the imag-
inary part. Note that the contributions of λ = 0 and 2 are
always negligible. See text for more details.
In Section IIA we used an approximation of negligi-
ble dependence of MDFF on the j, l indices (see Eq.20).
Generally, as the wave vector k(j,l) for each Bloch wave
changes slightly by an amount given by the corresponding
eigenvalue γ(j,l), the values of qz and q
z would change ac-
cordingly and therefore we should not be allowed to take
MDFF out of the sum over the indices j, l in the Eq. (24).
However, the change in qz (and q
z) induced by the eigen-
values γ(j,l) is small and can be neglected with respect
to the qz = χ0E/2E0 given by the energy loss E
26. To
demonstrate this we plot the dependence of MDFF on
qz, q
z for qx and qy corresponding to the main DFF and
MDFF terms, see Fig. 1. If qz is given in a.u.
−1 (atomic
units, 1 a.u.= 0.529178Å), typical values for L2,3 edges
of Fe, Co and Ni are around tenth of a.u.−1, whereas
typical values of γ(j,l) for strongly excited Bloch waves
are one or two orders of magnitude smaller. Thus the
approximation of weak j, l dependence of MDFF is well
justified.
Besides γ(j,l), the other factors determining the value
of qz are the energy of the edge, i.e. the energy lost by
the probe electron, the tilt with respect to the zone axis
and whether the excited beam is in a HOLZ. These last
factors have been included in our calculation. Only the
variations due to γ(j,l) are neglected, thus giving rise to
an error . 1%. If a more accurate treatment would be
needed, the smooth behavior of MDFF with respect to
qz would allow to use simple linear or quadratic interpo-
lation/extrapolation methods.
As mentioned in the introduction and explained in
Refs. 1,7, dichroism in the TEM is made possible by
the analogous role that the polarization vector ε and the
wave vector transfer q play in the dipole approximation
of the DDSCS. However we do not restrict our calcu-
lations to the dipole approximation. We use the more
complete expression Eq. (5).
To evaluate the accuracy of the dipole approximation,
we compare the dipole approximation of MDFF with the
full calculation (with λ up to 3) also showing λ-diagonal
components of the MDFF, Fig. 2. Because the domi-
nant contribution to the signal originates from (dipole
allowed) 2p→ 3d transitions, the λ = λ′ = 1 term nearly
coincides with the total MDFF. While the dipole approx-
imation works relatively well for the studied systems,
particularly the MDFF divided by squares of momen-
tum transfer vectors (right column of the Fig. 2), it has
significantly different asymptotic behaviours for larger q-
vectors. The λ = λ′ = 1 term provides a much better
approximation, which remains very accurate also in the
large q region.
It is worth mentioning that thanks to the properties
of the Gaunt coefficients the 2p → 3d transitions are all
included in the λ = 1 and λ = 3 contributions. Thanks
to the negligible value of the radial integrals for λ =
3 the terms with λ = 1 account for the large majority
of the calculated signal. The contributions from λ =
0, 2 describe transitions from 2p to valence p or f states
and are always negligible due to the composition of the
density of states beyond the Fermi level. They practically
overlap with the zero axis in all the six parts of Fig. 2.
It can be shown24 that in the dipole approximation
the real part of the MDFF is proportional to q · q′ and
the imaginary part is proportional to q × q′. A little
algebra can thus show that the imaginary part of the
MDFF is, in the geometry described in the caption of
Fig. 2, constant with respect to qx. As expected, the DFF
(which is proportional to q2) has a minimum at qx = 0,
where S(q, E)/q4 has a maximum. For the MDFF (and
corresponding S(q,q′, E)/q2q′2) these minima and max-
ima are centered at qx = −G/2 = −0.76 a.u.
−1 where
|qx| = |q
III. RESULTS
We summarize the results obtained for body-centered
cubic iron (bcc-Fe), hexagonal close-packed cobalt (hcp-
Co) and face-centered cubic nickel (fcc-Ni) crystals,
which are also the first samples prepared for EMCD mea-
surements. These results are valuable for optimization of
the experimental setup.
The geometry setup for observing the dichroic effect7
consists in creating a two-beam case by tilting the beam
away from a zone axis (here (001)) by a few degrees and
then setting the Laue circle center equal to G/2 for the
G vector to be excited. In analogy to XMCD, where two
measurements are performed for left- and righ-handed
circularly polarized light, here we perform two measure-
ments by changing the position of the detector, which
lies once at the top and once at the bottom of the Thales
circle having as diameter the line connecting the diffrac-
tion spots 0 and G. This geometry setup, together with
the crystal structure, is an input for the calculation of
the Bloch wave coefficients (within the systematic row
approximation) using the dynamical diffraction theory
code described in section IIA.
The electronic structure was calculated using the
WIEN2k package15, which is a state-of-the-art implemen-
tation of the full-potential linearized augmented plane
waves method. The experimental values of lattice param-
eters were used. More than 10000 k-points were used to
achieve a very good converge of the Brillouin zone inte-
grations. Atomic sphere sizes were 2.2, 2.3 and 2.2 bohr
radii for bcc-Fe, hcp-Co and fcc-Ni, respectively. The re-
sulting electronic structure was the input for the calcula-
tion of the individual MDFFs required for the summation
(see Section II B).
In the three studied cases the dichroic effect is domi-
nated by the transitions to the unoccupied 3d states. The
d-resolved spin-up density of states (DOS) is almost fully
occupied, while the spin-down d-DOS is partially unoccu-
pied. In Fig. 3 we compare the d-DOS with the dichroic
signal at the L3 edge. Due to negligible orbital moments
in these compounds the L2 edge shows a dichroic signal
of practically the same magnitude but with opposite sign.
The shape of the calculated dichroic peaks corresponds to
the difference of spin-up and spin-down d-DOS, similarly
to XMCD, as it was shown for the same set of systems in
Ref. 25. The calculations were performed within system-
atic row conditions with G = (200) for bcc-Fe and fcc-Ni
and G = (100) for hcp-Co. The sample thicknesses were
set to 20 nm, 10 nm and 8 nm for bcc-Fe, hcp-Co and fcc-
Ni, respectively. These values were found to be optimal
for these systems in the given experimental geometry.
An interesting point is the comparison of the strength
of the dichroic signal. According to the d-DOS projec-
tions one would expect comparable strength of signals
for the three elements under study. But the dichroic
signal of hcp-Co seems to be approximately a factor of
two smaller than that of the other two. The reason for
that can be explained by simple geometrical considera-
-6 -4 -2 0 2 4
E - E
 ( eV )
0 2 4
 ( eV )
FIG. 3: (online color) Spin-resolved d-densities of states (left) and resulting signal on L3 edge (right) on bcc-Fe, hcp-Co and
fcc-Ni (from top to bottom) at optimal thickness (see text). Spin-up DOS is drawn using a solid black line (positive) and spin-
down DOS using a dashed red line (negative). DDSCS for the (+) detector position is drawn using a solid blue line, DDSCS
for the (-) position is drawn using a dashed green line. The dichroic signal (difference) is the hatched red area. G = (200) for
bcc-Fe and fcc-Ni and (100) for hcp-Co.
tions starting from Eq. (24). For simplicity we consider
only the main contributions: the DFF S(q,q, E) and
the MDFF S(q,q′, E) with q ⊥ q′. For bcc-Fe and fcc-
Ni the summation over u within the Bravais cell leads
always to the structure factor 2 and 4, respectively, be-
cause q′ − q = G is a kinematically allowed reflection.
This factor cancels out after division by the number of
atoms in the Bravais cell. Therefore it does not matter,
what is the value of q-vectors, the sum over the atoms
is equal to S(q,q′, E)/q2q′2 itself. On the other hand,
the unit cell of hcp-Co contains two equivalent atoms at
positions u1 = (
) and u2 = (
). For the two
DFFs q = q′ and the exponential reduces to 1; since there
are two such terms, after division by Nu the sum equals
again the DFF itself. But for the main MDFF we have
q ⊥ q′ and the exponential factor will in general weight
the terms. One can easily see, that q′ − q = G. For the
G = (100) systematic row case, which was used for cal-
culation of hcp-Co in Fig. 3 the exponentials evaluate to
the complex numbers − 1
and − 1
for u1 and
u2, respectively. Because of symmetry, the MDFFs for
both atoms are equal and then the sum 1
leads to a
factor − 1
for the MDFF contribution, i.e. the influence
of its imaginary part, which is responsible for dichroism,
on the DDSCS is reduced by a factor of two.
To optimize the dichroic signal strength of hcp-Co, we
require G · u1 = G · u2 = 2πn, which gives in principle
an infinite set of possible G vectors. The one with lowest
hkl indices is G = (110). A calculation for this geometry
setup leads to approximately twice the dichroic signal,
see Fig. 4 and compare to the corresponding graph in
Fig. 3.
For the optimization of the experimental setup it is
important to know how sensitive the results are to varia-
tion of the parameters like the thickness of the sample or
the accuracy of the detector position. Another question
related to this is also the sensitivity to the finite size of
the convergence and collection angles α and β. In the
following text we will address these questions.
The thickness influences the factor Tjlj′l′ in the
Eq. (24) only. This factor leads to the so called pen-
dellösung oscillations - modulations of the signal strength
-1 0 1 2 3 4
Energy loss ( eV )
dichroism (%)
DDSCS(+)
DDSCS(-)
FIG. 4: L3 peak of hcp-Co calculated for the G = (110)
systematic row at 18nm. See caption of Fig. 3. The peaks
have been renormalized so that their sum is 1, therefore their
difference is the dichroic signal (ca. 15% in this case).
as a function of thickness. This also influences the
strength of the dichroic signal. Results of such calcula-
tions are displayed in Fig. 5 (we did not include absorp-
tion into these simulations, so that all signal variations
are only due to the geometry of the sample). From these
simulations it follows that a well defined thickness of the
sample is a very important factor. Relatively small vari-
ations of the thickness can induce large changes in the
dichroic signal, particularly in fcc-Ni. From the figure
one can deduce that the optimal thickness for a bcc-Fe
sample should be between 8 nm and 22 nm (of course, due
to absorption, thinner samples within this range would
have a stronger signal), for hcp-Co between 15 nm and
22 nm and for fcc-Ni it is a relatively narrow interval -
between 6 nm and 10 nm. However, we stress that these
results depend on the choice of the systematic row vector
G. For example hcp-Co with G = (100) (instead of (110)
shown in Fig. 5) has a maximum between 5 nm and 15
nm (although it is much lower, as discussed before).
Taking the optimal thickness, namely 20 nm, 18 nm
and 8 nm for bcc-Fe, hcp-Co and fcc-Ni, respectively, we
calculated the dependence of the dichroic signal on the
detector position. We particularly tested changes of the
dichroic signal when the detector is moved away from its
default position in the direction perpendicular to G, see
Fig. 6. It is interesting to note that the maximum abso-
lute difference occurs for a value of qy smaller than |G|/2.
This can be qualitatively explained by considering the
non-zero value of qz and q
z, i.e. q and q
′ are not exactly
perpendicular at the default detector positions. More-
over the MDFF enters the summation always divided by
q2q′2 and the lengths of q-vectors decrease with decreas-
ing qy. The important message we can deduce from this
figure is that the dichroic signal is only weakly sensitive
to the accuracy of qy since even displacement by 10-20%
in the detector default qy positions (i.e. qy = ±G/2) do
not affect significantly the measured dichroic signal.
0 20 40 60 80 100
Thickness ( nm )
FIG. 5: (online color) Dependence of the DDSCS and of the
dichroic signal on sample thickness for a) bcc-Fe, b) hcp-Co
and c) fcc-Ni. Systematic row vector G = (200) was used for
bcc-Fe and fcc-Ni, while for hcp-Co G = (110) was chosen.
The blue and green solid curves are DDSCSs calculated for
the (+) and (-) detector positions, the dashed black curve is
the DFF part of the DDSCS (it is identical for both detector
positions). The red line with circles is the relative dichroism
defined as difference of DDSCSs divided by their sum, the red
solid curve is the absolute dichroism - difference of DDSCSs.
Related to this is a study of the dependence of the
dichroic signal on the finite size of the convergence and
collection angles α and β. We performed a calculation
for the three studied metals and found that collection
and convergence half-angles up to 2 mrad weakens the
relative dichroic signal by less than 10%.
IV. CONCLUSIONS
We have developed a computer code package for the
calculation of electron energy loss near edge spectra,
which includes the theory of dynamical Bragg diffrac-
tion. We applied the code to the recently discovered phe-
nomenon of magnetic chiral dichroism in the TEM and
we demonstrated the relation of the dichroic peak shape
to the difference of d-projections of the spin-resolved
density of states in analogy with similar observation for
XMCD.
Using this code we examined the validity of the dipole
approximation, which is often assumed. We found that
-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8
 / |G| ( dimensionless )
bcc-Fe
hcp-Co
fcc-Ni
absolute dichroism (arb. units)
FIG. 6: (online color) Dependence of the dichroic signal on
detector displacements along qy . The full symbols correspond
to the relative dichroic signal while the open symbols to the
difference of the DDSCS for both detector positions. These
are in arbitrary units and their magnitudes are not directly
comparable. Vertical lines are showing the default detector
positions.
for the 3d ferromagnetic systems studied it is a reasonable
approximation, however with wrong asymptotic proper-
ties - it overestimates the contributions from larger q-
vectors. A very accurate approximation for the studied
systems is the λ = λ′ = 1 approximation, which treats
appropriately the dominant p→ d dipole transitions and
remains very accurate also for large q, q′.
In order to provide guidance to the experimentalist we
have investigated the strength of the dichroic signal as a
function of the sample thickness and the precision of the
detector placement. While the dichroic signal strength is
rather robust with respect to the precision of the detector
placement, the thickness of the specimen influences the
signal considerably. Therefore it might be a challenge to
produce samples with optimum thickness and selecting
the best systematic row Bragg spot. Our calculations
yield best thicknesses in order to detect EMCD of the
iron and nickel samples for the systematic rowG = (200)
to be 8-22 nm and 6-10 nm, respectively, and for cobalt
in the systematic row G = (110) to be 15-22 nm.
Acknowledgments
We thank Dr. Cécile Hébert and Dr. Pavel Novák for
stimulating discussions. This work has been supported
by the European Commission, contract nr. 508971 (CHI-
RALTEM).
∗ Electronic address: rusz@mag.mff.cuni.cz
† currently at Department of Physics, Uppsala University,
Box 530, S-751 21 Uppsala, Sweden
1 A. P. Hitchcock, Jpn. J. Appl. Phys. 32(2), 176 (1993).
2 J. Yuan and N. K. Menon, J. Appl. Phys. 81(8), 5087
(1997).
3 B. T. Thole, P. Carra, F. Sette and G. van der Laan, Phys.
Rev. Lett. 68, 1943 (1992).
4 S. W. Lovesey and S. P. Collins, X-Ray Scattering and Ab-
sorption by Magnetic Materials, Clarendon Press, Oxford,
UK, 1996.
5 J. Stöhr, J. Electron Spectrosc. Relat. Phenom. 75, 253
(1995).
6 C. Hébert and P. Schattschneider, Ultramicroscopy 96, 463
(2003).
7 P. Schattschneider, S. Rubino, C. Hébert, J. Rusz, J.
Kuneš, P. Novák, E. Carlino, M. Fabrizioli, G. Panaccione
and G. Rossi, Nature 441, 486 (2006).
8 M. Nelhiebel, Ph.D. Thesis, 1999, Vienna University of
Technology.
9 M. Nelhiebel, P. Schattschneider and B. Jouffrey, Phys.
Rev. Lett. 85(9), 1847 (2000).
10 P. Schattschneider, C. Hébert, H. Franco and B. Jouffrey,
Phys. Rev. B 72, 045142 (2005).
11 H. Kohl and H. Rose, Advances in electronics and electron
optics 65, 173 (1985).
12 P. Schattschneider, M. Nelhiebel and B. Jouffrey, Phys.
Rev. B 59, 10959 (1999).
13 L.-M. Peng and M. J. Whelan, Proc. R. Soc. London Ser.
A 431, 111 (1990).
14 A. J. F. Metherell, in Electron Microscopy in Materials
Science II, 397 (1975), CEC, Luxembourg.
15 P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka and
J. Luitz, 2001 WIEN2k, Vienna University of Technology
(ISBN 3-9501031-1-2).
16 The WIEN2k package can calculate X-ray structure fac-
tors. By supplying the potential instead of charge density
it is possible to use the same code to calculate the electron
structure factors using a simple script.
17 A. Weickenmeier and H. Kohl, Acta Cryst. A47, 590
(1991).
18 P. A. Doyle and P. S. Turner, Acta Cryst. A24, 390 (1968).
19 A. L. Lewis, R. E. Villagrana and A. J. F. Metherell, Acta
Cryst. A34, 138 (1978).
20 Y. Kainuma, Acta Cryst. 8, 247 (1955).
21 P. Schattschneider and W. S. M. Werner, J. Electron Spec-
trosc. Relat. Phenom. 143, 81 (2005).
22 P. Hohenberg, W. Kohn, Phys. Rev. 136, B864 (1964) and
W. Kohn, L. J. Sham, Phys. Rev. 140, A1133 (1965).
23 J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
24 P. Schattschneider, C. Hébert, S. Rubino, M. Stöger-
Pollach, J. Rusz, P. Novák, submitted.
25 J. Kuneš and P. M. Oppeneer, Phys. Rev. B 67, 024431
(2003).
26 The question of momentum conservation in the z direction
in the inelastic interaction in a crystal of finite thickness is
related to the probability of inter- and intrabranch transi-
tions of the probe electron?
mailto:rusz@mag.mff.cuni.cz
