{'fraction_non_alphanumeric': 0.03305640270620815, 'fraction_numerical': 0.0116030427314474, 'mean_word_length': 4.302602230483271, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
First principles theory of chiral dichroism in electron microscopy applied to 3d ferromagnets

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.


I. INTRODUCTION
The analogy between X-ray absorption spectroscopy (XAS) and electron energy loss spectroscopy (EELS) has been recognized long ago 1,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 counterpart of X-ray magnetic circular dichroism (XMCD) 3,4,5 experiments with electron probes was thought to be technically 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 electron with l = 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 circular dichroism using a standard non-polarized electron beam in the TEM 6 it was demonstrated that such experiments (called energy-loss magnetic chiral dichroism, EMCD) are indeed possible 7 . This novel technique is of considerable interest for nanomagnetism and spintronics 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 approximations for the mixed dynamic form factor. In section III we study the dependence of the dichroic signal of bcc-Fe, hcp-Co and fcc-Ni on various experimental conditions. 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 approximation 10 DDSCS is written as with where q = k f −k 0 is the difference (wave vector transfer) between final wave vector k f and initial wave vector k 0 of the fast electron; γ = 1/ 1 − v 2 /c 2 is a relativistic factor and a 0 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 superposition of Bloch waves, which reflects the discrete translation symmetry. Each Bloch wave can be decomposed into a linear combination of plane waves -it is a coherent superposition of (an in principle infinite number of) plane waves. The wave function of the fast electron can be thus written as for incident wave and for outgoing wave, where C h are so called Bloch coefficients, ǫ (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 generalization of the DFF -the mixed dynamic form factors 11 (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 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 jl gh = k (l) − k (j) + h − g and the total DDSCS will be a sum over all diads of q and q ′ vectors of terms 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) coefficients are given by dynamical diffraction theory. This will be covered in the next subsection II A. 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 dynamical diffraction theory and identification of important terms. This task is mainly 'geometry dependent', although it can also contain some input from electronic structure codes, namely the Coulomb part of crystal potential. ii) Calculation of MDFFs requested by the dynamical 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 generalization of the formalism presented in Ref. 8,12 extending 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 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 Schrödinger equation we obtain the secular equation 14 where K 2 = U 0 + 2meE/ 2 , m and e are, respectively, the electron mass and charge, U g = 2meV g / 2 where V g are the Fourier components of the crystal potential, which can be either calculated ab initio 15,16 or obtained from the tabulated forms of the potential 17,18 . It can be shown 14,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 matrix 13,19 This eigenvalue problem can be transformed into a hermitean one using a diagonal matrix D with elements Then the eigenvalue problem is equivalent to and the original Bloch wave coefficients can be retrieved using the relation By solving this eigenvalue problem we obtain the fast electron wave function as a linear combination of eigenfunctions 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 crystal surface. The crystal surface is a plane defined by the scalar product n · r = t 0 . Then the boundary condition (in the high energy limit) leads to the following condition 14 It is easy to verify that as required by the boundary condition. We have used the continuity condition, Eq. (7), and the completeness relation for the Bloch coefficients Therefore the wave function of the fast moving electron in the crystal, which becomes a single plane wave at n·r = t 0 is given by the following expression The following discussion will be restricted to a particular case -a crystal with parallel surfaces. For such a crystal with normals in the direction of the z axis we set t 0 = 0 for the fast electron entering the crystal and t 0 = 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 electron wave function (Bloch field) onto a plane wave after the inelastic event. Therefore the calculation of the ELNES requires the solution of two independent eigenvalue problems describing an electron wave function before and after the inelastic event 8,12 . Invoking reciprocity for electron propagation the outgoing wave can also be considered as a time reversed solution of the Schrödinger equation, also known as the reciprocal wave 20 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 (j) g for the Bloch coefficients of the incoming electron and we use D (l) h for the Bloch coefficients of the outgoing electron entering the detector (obtained from the two independent eigenvalue problems). Similarly, superscript indices (j) and (l) indicate eigenvalues and Bloch-vectors for incoming and outgoing electron, respectively. We thus obtain In crystals the position of each atom can be decomposed 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 approximation 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 vectors. The sum in Eq. 6 can be separated into two terms Since and the algebric sum of g, h is simply a reciprocal lattice vectors G, which fulfills e iG·R = 1, it is possible to simplify the second term 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 leads 8,21 to so that the total sum over all atomic positions is ). The final expression of the DDSCS we write as where depends only on the eigenvectors of the incoming and outgoing beam and 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 outgoing wave, respectively. Particular η (j) can be calculated using the following expression 14 and similarly for the outgoing beam. This way the eigenvalues change from γ (j) to γ (j) +iη (j) and the ∆ acquires an imaginary part. Such approximative treatment of absorption thus affects only the thickness-dependent function T jlj ′ 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 10 8 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 considerable 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 U g . A dimensionless parameter w g = s g ξ g -product of the excitation error and the extinction distance 14 -reflects both these criteria. Therefore we can filter the list of beams by selecting only beams with w g < w max . Experience shows that in the final summation 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 (several hundreds). Therefore we defined two cut-off parameters for w g -the first for the solution of the eigenvalue problem (typically w max,1 is between 1000 and 5000) and the second for the summation (w max,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 product of their excitation ǫ (j) and their norm on the subspace defined by selected subset of beams, C (j) 0 ||C (j) || subsp . In the systematic row conditions this value is large only for a small number of Bloch waves. Typically in the experimental 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 between 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, orbital 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 classically, which leads to somewhat different expression for MDFF giving incorrect L 2 − L 3 branching ratio. The final expression is 7 Here we made use of Wigner 3j-symbols, Y λ µ are spherical 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 D LMS (νk) is the projection of the (νk) Bloch state onto the LM S 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 D LMS (νk) we employ the density functional theory 22 within the local spin density approximation 23 . In Section II A we used an approximation of negligible 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 q z and q ′ z would change accordingly 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 q z (and q ′ z ) induced by the eigenvalues γ (j,l) is small and can be neglected with respect to the q z = χ 0 E/2E 0 given by the energy loss E 26 . To demonstrate this we plot the dependence of MDFF on q z , q ′ z for q x and q y corresponding to the main DFF and MDFF terms, see Fig. 1. If q z is given in a.u. −1 (atomic units, 1 a.u.= 0.529178Å), typical values for L 2,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 q z 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 q z would allow to use simple linear or quadratic interpolation/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 calculations 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 dominant contribution to the signal originates from (dipole allowed) 2p → 3d transitions, the λ = λ ′ = 1 term nearly coincides with the total MDFF. While the dipole approximation works relatively well for the studied systems, particularly the MDFF divided by squares of momentum transfer vectors (right column of the Fig. 2), it has significantly different asymptotic behaviours for larger qvectors. 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 shown 24 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 q x . As expected, the DFF (which is proportional to q 2 ) has a minimum at q x = 0, where S(q, E)/q 4 has a maximum. For the MDFF (and corresponding S(q, q ′ , E)/q 2 q ′2 ) these minima and maxima are centered at q x = −G/2 = −0.76 a.u. −1 where |q x | = |q ′ x |.

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 measurements. These results are valuable for optimization of the experimental setup.
The geometry setup for observing the dichroic effect 7 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 measurements 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 diffraction 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 II A.
The electronic structure was calculated using the WIEN2k package 15 , which is a state-of-the-art implementation of the full-potential linearized augmented plane waves method. The experimental values of lattice parameters were used. More than 10000 k-points were used to achieve a very good converge of the Brillouin zone integrations. Atomic sphere sizes were 2.2, 2.3 and 2.2 bohr radii for bcc-Fe, hcp-Co and fcc-Ni, respectively. The resulting electronic structure was the input for the calculation of the individual MDFFs required for the summation (see Section II B).
In the three studied cases the dichroic effect is dominated 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 unoccupied. In Fig. 3 we compare the d-DOS with the dichroic signal at the L 3 edge. Due to negligible orbital moments in these compounds the L 2 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 systematic 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 projections 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- 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, because 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)/q 2 q ′2 itself. On the other hand, the unit cell of hcp-Co contains two equivalent atoms at positions u 1 = ( 1 3 , 2 3 , 1 4 ) and u 2 = ( 2 3 , 1 3 , 3 4 ). For the two DFFs q = q ′ and the exponential reduces to 1; since there are two such terms, after division by N u 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 calculation of hcp-Co in Fig. 3 the exponentials evaluate to the complex numbers − 1 2 ± i Nu u leads to a factor − 1 2 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 · u 1 = G · u 2 = 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 variation 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 T jlj ′ l ′ in the Eq. (24) only. This factor leads to the so called pendellösung oscillations -modulations of the signal strength as a function of thickness. This also influences the strength of the dichroic signal. Results of such calculations are displayed in Fig. 5 (we did not include absorption 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 variations 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 intervalbetween 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 absolute difference occurs for a value of q y smaller than |G|/2. This can be qualitatively explained by considering the non-zero value of q z and q ′ z , i.e. q and q ′ are not exactly perpendicular at the default detector positions. Moreover the MDFF enters the summation always divided by q 2 q ′2 and the lengths of q-vectors decrease with decreasing q y . The important message we can deduce from this figure is that the dichroic signal is only weakly sensitive to the accuracy of q y since even displacement by 10-20% in the detector default q y positions (i.e. q y = ±G/2) do not affect significantly the measured dichroic signal. 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 diffraction. We applied the code to the recently discovered phenomenon 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 for the 3d ferromagnetic systems studied it is a reasonable approximation, however with wrong asymptotic properties -it overestimates the contributions from larger qvectors. 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 row G = (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.