Citation |

- Permanent Link:
- http://ufdc.ufl.edu/AA00029942/00001
## Material Information- Title:
- Electron transport theory in magnetic nanostructures
- Creator:
- Choy, Tat-Sang, 1971-
- Publication Date:
- 2001
- Language:
- English
- Physical Description:
- viii, 123 leaves : ill. ; 29 cm.
## Subjects- Subjects / Keywords:
- Atoms ( jstor )
Conductivity ( jstor ) Electrical resistivity ( jstor ) Greens function ( jstor ) Impurities ( jstor ) Magnetism ( jstor ) Magnets ( jstor ) Matrices ( jstor ) Superlattices ( jstor ) Vertices ( jstor ) Dissertations, Academic -- Physics -- UF ( lcsh ) Electron transport ( lcsh ) Physics thesis, Ph. D ( lcsh ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Thesis:
- Thesis (Ph. D.)--University of Florida, 2001.
- Bibliography:
- Includes bibliographical references (leaves 118-122).
- General Note:
- Printout.
- General Note:
- Vita.
- Statement of Responsibility:
- by Tat-Sang Choy.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. Â§107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
- Resource Identifier:
- 027677960 ( ALEPH )
48188190 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

ELECTRON TRANSPORT THEORY IN MAGNETIC NANOSTRUCTURES BY TAT-SANG CHOY A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2001 ACKNOWLEDGMENTS I would like to thank Selman Hershfield for being my thesis advisor. He introduced me to the field of magnetic nanostructures and taught me Green's function techniques. I also learned from Selman how to work with experimentalists. He has given me guidance in research and also encouraged me to find research topics myself. When we started working on optimizing the GMR structure, I was not optimistic that the difficulties could be overcome. Without Selman's encouragement and support, I would have waited a few years before taking this challenge. I would like to thank Professors Arthur Hebard, Peter Herschfeld, Mark Meisel, Stephen Pearton, Fred Sharifi, and Christopher Stanton for serving as as members of my supervisory committee. Fred Sharifi taught me a lot about sample growth, information I used in my research. I also enjoyed a lot the ad hoc discussion group organized by Professors Arthur Hebard, Selman Hershfield, Mark Meisel, Andrew Rinzler, and Fred Sharifi. Two good friends influenced me greatly in my physics career. Hoi-Fung Chau was the first research physicist I have known. Ever since high school, I was lucky to have learned from him many new ideas in physics, mathematics, and computer sciences. Those ideas helped me realize my true passion is in physics rather than computer science. Jian Chen taught me a lot during and after he worked in Selman's group as a postdoc. I enjoyed dropping by his office chatting with him about physics and other things. I regard Jian Chen as my first close collaborator in science. I would like to thank fellow graduate students Kingshuk Majumdar, Steven Patamia, Vladimir Boychev, Zhihong Chen, Stephanie Getty, Heather Hudspeth, ii 7 ' Quentin Hudspeth, and Kelvin McCarthy for their friendship. Steven Patamia gave me useful advice in C++ programming. Many thank to our former system administrator Chandra Chegireddy, who helped tremendously in setting up the special hardware and software for my research. I would like to thank Dr. Michael Jones for making the LaTex dissertation template available. The work in this thesis was supported by grants from the Air Force Office of Scientific Research, the National Science Foundation, and the Research Corporation Company. Finally, I would like to thank my parents and my sister for their support. I am grateful for my parents' understanding of their son spending the past few years eleven time zones away. Also, my sister's gifts from time to time are much appreciated. iii TABLE OF CONTENTS ACKNOWLEDGMENTS . . . . . . . . . . ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTERS 1 INTRODUCTION . . . . . . . . . . . . 1.1 1.2 1.3 1.4 GMR Experiments ................ Models of the CPP-GMR ............. 1.2.1 Resistors-in-Series Model .......... 1.2.2 Fermi Surface Mismatch Picture ..... 1.2.3 Potential Step Picture .......... 1.2.4 Effects of the Disorder .......... Optimization at the Atomic Level . . . . . . . . Overview of the Rest of the Dissertation . . . . . . . . . . . . . . . 2 . . . . . . . . . . . 4 . . . . . . . . . . . 5 . . . . . . . . . . . 7 . . . . . . . . . . . 10 . . . . . . . . . . . 12 . . . . . . . . . . . 14 . . . . . . . . . . . 17 2 MODELING TRANSPORT IN NANOSTRUCTURES . 2.1 Modeling of the Band Structure ............ 2.1.1 Multi-band Tight-binding Hamiltonian . ... 2.1.2 Modeling Impurities . . . . . . . . . . . . . . 2.1.3 Impurity Averaged Green's Functions . . . . . 2.2 Conductivity Calculation . . . . . . . . . . . . . . . . 2.2.1 Formal Aspects of Current Operators . . . . . 2.2.2 Current Operators in One Dimension . . . . . 2.2.3 Kubo Formula . . . . . . . . . . . . . . . . . . 2.2.4 Current Conservation . . . . . . . . . . . . . . 2.3 Realistic Parameters . . . . . . . . . . . . . . . . . . 2.3.1 Tight-Binding Parameters . . . . . . . . . . . 2.3.2 2.3.3 2.3.4 18 20 20 21 23 26 27 31 36 39 43 43 44 46 46 Impurity Potential for Substitutional Disorder Impurity Potential for Structural Disorder . . Extension to Magnetic Alloys . . . . . . . . . iv . . . . . . . . . . . ii vii . . 3 COMPUTATION PROCEDURE . . . . . . . . . . . . . . 3.1 Current-Perpendicular-to-Plane Conductivity ........ 3.1.1 Impurity Averaged Green's Function . . . . . . . . . 3.1.2 Vertex Correction . . . . . . . . . . . . . . . . . . . . 3.1.3 Conductivity . . . . . . . . . . . . . . . . . . . . . . 3.1.4 Diagonal Impurity Potential . . . . . . . . . . . . . . 3.1.5 Total Computational Cost . . . . . . . . . . . . . . . 3.2 Current-In-Plane Conductivity . . . . . . . . . . . . . . . . . 3.3 Conductivity Tensor . . . . . . . . . . . . . . . . . . . . . . 3.4 C oding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Template Numerical Toolkit (TNT) . . . . . . . . . . 3.4.2 3.4.3 Fortran performance with BLAS and LAPACK . . . BlockMatrix class . . . . . . . . . . . . . . . . . . 4 APPLICATION TO GIANT MAGNETORESISTANCE . . . . . . . 4.1 Experiment of Gijs and Collaborators . . . . . . . . 4.2 Experiment of Cyrille and Collaborators . . . . . . 4.3 Comparing Theory with Experiments . . . . . . . . 4.4 Other Effects on the GMR . . . . . . . . . . . . . . 5 OPTIMIZATION OF GMR STRUCTURES . . . . . . 5.1 5.2 5.3 5.4 5.5 5.6 5.7 . . . . . . . . . 70 . . . . . . . . . 71 . . . . . . . . . 74 . . . . . . . . . 81 84 Inverse Problems . . . . . . . . . . . . . . . . . . . . . . . . Simulated Annealing . . . . . . . . . . . . . . . . . . . . . . Accuracy and Speed Requirements on Model . . . . . . . . . Constrained Optimization of the GMR . . . . . . . . . . . . Optimal GMR Configuration in CoNiCu superlattices . . . . Understanding the Large GMR . . . . . . . . . . . . . . . . C onclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 . . . . 87 . . . . 88 . . . . 90 . . . . 92 . . . . 94 . . . . 97 6 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . 100 APPENDICES A FAST INVERSE ALGORITHMS . . . . . . . . . . . . . . 103 A.1 Godfrin Inverse Algorithm for Tridiagonal Block Matrix . . . . . . . 104 A.2 Inverse of Tridiagonal Block Matrix with Periodic Boundary Condition 108 v 49 . . . . 50 . . . . 51 . . . . 53 . . . . 54 . . . . 56 . . . . 58 . . . . 59 . . . . 60 . . . . 61 . . . . 62 . . . . 64 . . . . 64 69 A.2.1 Rank-one changes: Woodbury formula . . . . . . . . . . . . 108 A.2.2 Applying the Woodbury Formula . . . . . . . . . . . . . . . 109 A.2.3 Cost Considerations . . . . . . . . . . . . . . . . . . . . . . 110 A.3 Efficient Coding for the Godfrin Algorithm . . . . . . . . . . . . . . 112 B COHERENT POTENTIAL APPROXIMATION . . . . . . . . . 115 REFERENCES........ .......................... .. 118 BIOGRAPHICAL SKETCH . . . . . . . . . . . . . . . . . . . . 122 vi Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy ELECTRON TRANSPORT THEORY IN MAGNETIC NANOSTRUCTURES By Tat-Sang Choy August 2001 Chairman: Selman P. Hershfield Major Department: Physics Magnetic nanostructure has been a new trend because of its application in making magnetic sensors, magnetic memories, and magnetic reading heads in hard disks drives. Although a variety of nanostructures have been realized in experiments in recent years by innovative sample growth techniques, the theoretical study of these devices remain a challenge. On one hand, atomic scale modeling is often required for studying the magnetic nanostructures; on the other, these structures often have a dimension on the order of one micrometer, which makes the calculation numerically intensive. In this work, we have studied the electron transport theory in magnetic nanostructures, with special attention to the giant magnetoresistance (GMR) structure. We have developed a model that includes the details of the the band structure and disorder, both of which are both important in obtaining the conductivity. We have also developed an efficient algorithm to compute the conductivity in magnetic nanostructures. The model and the algorithm are general and can be applied to complicated structures. We have applied the theory to current-perpendicular-toplane GMR structures and the results agree with experiments. Finally, we have vii searched for the atomic configuration with the highest GMR using the simulated annealing algorithm. This method is computationally intensive because we have to compute the GMR for 103 to 104 configurations. However it is still very efficient because the number of steps it takes to find the maximum is much smaller than the number of all possible GMR structures. We found that ultra-thin NiCu superlattices have surprisingly large GMR even at the moderate disorder in experiments. This finding may be useful in improving the GMR technology. viii CHAPTER 1 INTRODUCTION Magnetic nanostructures [1-2] are magnetic devices which have at least one dimension of the order a nanometer, 10-9 m. These devices have applications in making magnetic sensors, magnetic memories, and magnetic reading heads in hard disks drives. Predicting the electrical transport properties theoretically is a challenge because the detailed arrangement of atoms affects the observable properties significantly. In this dissertation, we will develop methods to predict the transport properties of magnetic nanostructures by taking into account band structure and disorder. We will first develop models to compute the giant magnetoresistance ratio (GMR), which is the fractional change of the resistance due to an applied magnetic field, of multilayer structures and compare our calculations with experiments. We will then solve the "inverse problem," which is to find the atomic configuration that has the highest GMR. Although the method we develop can be applied to various geometries, in this dissertation we concentrate on magnetic superlattices, which are multilayers formed by repeatedly stacking magnetic layers followed by non-magnetic layers. The thickness of each layer is of the order of one nanometer; however, the other two dimensions can be as large as one micrometer, 10-6 m. This geometry can be grown by techniques used in depositing thin films. In the following, we will first discuss the GMR experiments, explain several pictures for understanding the gi- 1 2 ant magnetoresistance, and discuss the recent trend of optimization on the atomic level. 1.1 GMR Experiments The first "giant" magnetoresistance measurement was that of Baibich et al. [3], where they observe an almost 50% reduction in the resistance of Fe/Cr superlattices in an applied magnetic field of 20 kG. As shown in Figure 1.1, the resistance at zero field, RO, is the highest. As the magnetic field is applied, the magnetization of the adjacent Fe layers align, which causes the resistance to decrease. At a field larger than the saturation field Hs, the magnetization of the adjacent Fe layers is parallel to each other. Therefore, the resistance stays at the value of the saturation field, Rsat. The giant magnetoresistance ratio is defined as GMR = Ro - R(1t Rsat In this experiment, the current travels in the plane of the multilayers, which is referred to as the current-in-plane or CIP giant magnetoresistance. In this dissertation, we focus on another geometry, the current-perpendicularto-plane (CPP) geometry, in which the current flows perpendicular to the plane of the multilayer. The CPP-GMR is usually larger than the CIP-GMR for the same lattice. Most experiments to date are for the CIP geometry, while a smaller portion of the experiments are for the CPP geometry [4-15]. Current-perpendicular-toplane experiments are difficult to perform because the thickness of a superlattice is of the order of a hundred nanometers at most, while the area of the layer is 1 0 8 06 -40 -30 -20 -10 0s R/R(H=0) (Fe 30 A/ Cr 18 A)0 Hs (Fe 30 /Cr 12 A h5 Hs (Fe 30 A / Cr 9A)0 - Hs 0 10 20 30 10 Magnetic field (kG) Figure 1.1: Magnetoresistance of three Fe/Cr superlattice at 4.2 K [3]. The current and the applied field are along the same [110] axis in the plane of the layers. H, is the saturation field required to align the magnetizations of Fe layers. As the magnetizations are aligned, the resistance is reduced from the peak at H = 0. usually much larger than 1 pm x 1 pm. When a current passes though the system perpendicular to the plane, the resistance is too small compared with the resistance of the lead and the contact, and is thus hard to measure. To measure the CPP resistivity, one must improve the aspect ratio of the sample. Small cross-sectional area to length ratio has been produced by using micro-structured pillars (Figure 1.2) [4], by growing superlattices inside the pores in polycarbonate (Figure 1.3) [15], and by connecting several superlattices with superconducting leads (Figure 1.4) [6-7]. 4 Mo 100 nm Polyimide -,A1203 200 nm -Mo 100 nm - Au 300 nm Fe/Cr 500 nm Cr 400 nm substrate SiO2 (d) Au (h) AuFe/Cr Figure 1.2: Schematic diagram the samples growth process used by Gijs et al. [4]. The cross sectional area S ranges from 6 to 130 pm2 1.2 Models of the CPP-GMR There has been important progress in understanding the CPP-GMR effect. Classical resistors-in-series models can be applied to study the CPP-GMR when the layer thickness is much larger than the Fermi wavelength [16]. The origin of the spin anisotropy, which measures the difference between the resistivity in the two spins, is more difficult to understand. Various free-electron models [17-21], tightbinding models [22-26], and density-function methods [27] have been performed. In the following, we discuss several models for the CPP-GMR. In the CPPGMR system, each electron has to go through every layer and interface. It is thus easier to understand intuitively than CIP-GMR system, where electrons travels along the layers, while at the same time going across interfaces of layers and/or scattering from the interfaces. 5 40r m Polycarbonate Cu layer Co layer CU film Figure 1.3: Schematic representation of the array of nanowires in the insulating polymer matrix [15]. 1.2.1 Resistors-in-Series Model When the layer thickness is smaller than the spin-diffusion length [28-29], one can treat the two spins independently, and this approach is called the two-fluid model [16]. When the layer thickness is much larger than the Fermi wavelength, the GMR effect can be explained with the classical resistors-in-series model, which assumes electrons pass through spin-dependent resistors, as shown in Figure 1.5. Let the resistance for the majority and minority spin in the bilayer formed by a magnetic layer and a non-magnetic layer be R and r, respectively. At saturation field, the magnetic moment of adjacent magnetic layers are parallel to each other. Assuming the spin diffusion length is much larger than the layer thickness, the resistance for the majority and minority spin is R + R and r + r, as shown in 6 e Cr Nb SiO (a) (b) Figure 1.4: CPP sample of Cyrille et al. [6]. (a) Schematic cross section of a CPP sample; (b) optical micrograph of a typical sample. Figure 1.5(a). Therefore, the net resistance for the saturation field is given by 1lRr Rrat - r (1.2) 2 R+ r At zero field, the magnetic moment of adjacent magnetic layers are assumed to be anti-parallel to each other, as shown in Figure 1.5(b). A majority spin electron of the first magnetic layer will be in the minority spin channel of the second magnetic layer. Therefore, the resistance for each spin channel is r + R. As a result, the net 7 resistance at zero field is given by Ro = R+r (1.3) 2 Since (R+ r)2 > 4Rr, one can see that Ro is larger than Rsat, which is the explains the GMR. (a) (b) R r r R R r R r Figure 1.5: Simple resistor model for CPP-GMR. The circuits corresponding to (a) parallel, and (b) anti-parallel, magnetic configuration. In both cases, the majority and minority spin electrons pass through the left and right series of resistors, respectively. In (a), the resistance for the majority and minority spin channels are 2R and 2r, respectively; the net conductance is 1/2R+ 1/2r. In (b), the resistance for either channels are r + R; the net conductance is 2/(r + R). 1.2.2 Fermi Surface Mismatch Picture Although the basic idea of the GMR can be understood by the simple resistor model above, this model does not explain the origin of the spin anisotropy, which is the difference divided by the sum of the resistance in the two spin channels. There are two spin anisotropies, the bulk spin anisotropy and the interface spin anisotropy. The interface spin anisotropy can be understood with simple pictures 6 and is often the dominant effect of the GMR unless the layers are thick. Therefore, in the following we will explain the GMR using only the interface spin anisotropy. The most intuitive way of explaining the spin anisotropies and the CPP-GMR is probably the Fermi surface picture. We use the Fermi surfaces plots from our "Fermi Surface Database" website at http: //www. uf 1. edu/fermisurf ace [3031] in the following discussion. It provides both pictures and interactive Virtual Reality Modeling Language (VRML) models of Fermi surfaces for most elemental metals. We assume the layers to be large compared with the Fermi wavelength such that the bulk Fermi surfaces is a reasonable approximation. Again the two fluid model is used, assuming the spin diffusion length is large compared to the thickness of the layers [16, 28-29]. As an example, we consider the Fe/Cr/Fe system, a three-layered structure formed by sandwiching a Cr layer in between two Fe layers. The Fermi surfaces of bulk Fe spin-up and spin-down channels, and non-magnetic Cr, are arranged in Figure 1.6 for both (a) the parallel and (b) the anti-parallel magnetic configurations. As shown in Figure 1.6(a), the Cr Fermi surface is very different from the Fe spin-up Fermi surface, while it is much closer to the Fe spin-down Fermi surface. If an electron in the spin-up channel of the Fe layer on the left needs to go into the Cr layer, there will be larger scattering at the interface because a wavevector on the spin-up Fermi surface of bulk Fe is not on the bulk Cr Fermi surface. When this electron goes from the Cr layer to the Fe layer on the right, it is scattered again. Therefore, an electron goes though two large scattering events, and as a result, the spin-up channel has a low conductivity. On the other hand, a spin-down electron from the Fe layer on the left is not scattered as much in going 9 through the two interfaces. The scattering is much smaller because the Cr and Fe spin-down Fermi surfaces are similar. Thus, the spin-down channel has a higher conductivity. The total conductivity for the parallel configuration is the sum of the two channels, which is dominated by the spin-down channel. (a) (b) F , Cr FeupSpin Up Channel Fedn Cr edn Spin Down Channel. Fe Cr Fedn Spin Up Channel Fedn Cr Feu Spin Down Channel Figure 1.6: Fermi Surface Mismatch in the Fe/Cr/Fe multilayer. The bulk Fermi for each spin in the layers in the (a) parallel, and (b) anti-parallel, magnetic configuration are shown. There are only three different Fermi surfaces, Fe majority spin, Fe minority spin, and Cr. Since the spin of an electron crossing the layers is assumed to fixed, the electron behave as a majority or minority electron depending on whether its spin is aligned or anti-aligned with the local magnetization of the layer. The anti-parallel case is shown in Figure 1.6(b). The spin-up and spin-down Fermi surfaces for the Fe layer on the right are exchanged. When the local magnetization on the left and the right are anti-parallel to each other, an electron which 10 starts as a spin-up (majority) electron on the left will be in the minority population on the right. Hence it behave as a spin-down electron in the bulk of the Fe layer on the right. Here, electrons in the both spin channel go though the same amount of scattering. Since the Cr/Fe spin-up Fermi surface mismatch is much larger than the Cr/Fe spin-down Fermi surface mismatch, the conductivity for the spin-down channel of the parallel configuration is much larger than the other channels. Therefore, the conductivity for the parallel magnetic configuration is much larger than the conductivity for the anti-parallel magnetic configuration. In this case, the giant magnetoresistance is positive. The CPP-GMR in Co/Cu/Co can be explained similarly using the Fermi surfaces in Figure 1.7. Here, the Cu Fermi surface is similar to the Co spin-up Fermi surface, but rather different from the Co spin-down Fermi surface. Therefore, the conductivity is dominated by the spin-up channel of the parallel configuration. Again, the conductivity for the parallel magnetic configuration is significantly larger than the conductivity for the anti-parallel magnetic configuration. Hence, the GMR is also positive. We note that the spin anisotropy for the Fe/Cr interface has a different sign than the spin anisotropy for the Co/Cu interface. 1.2.3 Potential Step Picture A related way of understanding the spin anisotropy, or the difference in resistance for difference spins, at interfaces is to look at the potential step across the interfaces. As will be explained in Chapter 2, the major difference in the 11 (a) (b) co- tcup Cu Co-fccup Spin Up Channel ofL Cu Cofj~ Spin Down Channel O fccup CU cn Spin Up Channel CoSpnDwnCa Cu COfccup ISpin Down Channel Figure 1.7: Fermi Surface Mismatching the Co/Cu/Co multilayer. The bulk Fermi for each spin in the layers in the (a) parallel, and (b) anti-parallel, magnetic configuration are shown. tight-binding parameters in the 3d transition metals is the on-site energy of the d-orbitals. The system looks like a simple transmission-though-potential-barrier problem. We must be careful, however, that the potential step affects only the d-electrons, while the s- and p-electrons carry a significant fraction of the current. The effect is complicated by the spd-hybridization. From the tight-binding parameters, the Cu potential is 0.82 eV lower than the spin-up Co potential and is 2.3 eV higher than spin-down Co potential. Since the potential step in Cu/Co spin-up interface is higher, it has higher resistance. 12 Similarly, the Cr potential is 1.53 eV higher than the spin-up Fe potential and is 0.23 eV lower than the spin-down Fe potential. The potential step may also be used to estimate the impurity potential due to interface diffusion of, say Cu into Co. 1.2.4 Effects of the Disorder Although the simple pictures in the last section can help understand the origin of the giant magnetoresistance effect, these pictures do not include both the band structure effect and the disorder effect. We have explained that the band structure is often the source of the GMR. On the other hand, both the conductivity and the size of the GMR are determined by the disorder. If we compute the conductivity with small disorder, we would get a conductivity orders of magnitudes larger than what is usually obtained in experiments. Also, if the disorder is too small, the GMR would be too large compared to experiments. Therefore, disorder must be included in the calculation in order to obtain the conductivity and the GMR. As explained by J. Chen et al. [25], the disorder can affect the the density of states (DOS) significantly. The density of states of Cu and Co with different disorder parameters is shown in Figure 1.8. As we can see, the disorder significantly changes the density of states. The effect of the disorder is to spread out the DOS. The net effect at the Fermi surface is hard to predict without the calculations done in Chapter 2. In Cu and the spin-down channel of Co, the DOS at the Fermi level is not significantly affected. However, for the spin-up channel of Co, which has a similar band structure to Cu, the DOS is doubled due to the DOS spread from 13 the d-peak just below the Fermi level. We will see in the next chapter that this significantly reduces the conductivity, even though the DOS is doubled. 2.0 - n,=0.01 1.5 - ,- ---r.,=0.08 1.0 Cu 0.5 Eo.o 22.0 > 1.5 a1.0 - Cot fa-0.5 C/) 0.0 I I 2.0 Col 1.5 1.0 0 0.5% -4.0 -2.0 0.0 2.0 E-EF W) Figure 1.8: Density of states (DOS) of Cu, and Co spin-up and spin-down channels for different impurity density. When the impurity density is small, the DOS remain close to the unperturbed DOS. However, when the impurity density is moderate, the DOS can change a lot. The conductivity is affected by the DOS near the Fermi level. For Cu and Co spin-down, the change in DOS at the Fermi level are small. However, for Co spin-up channels, which has a similar band structure to Cu, the DOS at the Fermi level is increased by a factor of two [25]. To calculate the GMR, we have to find the appropriate disorder parameters. We set the bulk impurity parameter so that the bulk resistivity is similar to experiments. A single bulk impurity parameter should be sufficient to give the right bulk resistivity for all elements in a sample. In Chapter 1, we explained how to find 14 the impurity potential due to diffusion. The density of the impurity is harder to set because there is no experiment the clearly measure this parameter. However, since only the total scattering rate matters, one can set the parameters to match the experiments. Since the same impurity density has to be used for both parallel and anti-parallel magnetic configurations, the single parameter should be able to give the correct resistivity for both configuration. In other words, there should less parameter than the number of data points. In Figure 1.9, the current-perpendicular-to-plane resistivity times the layer thickness t for the Cot/Cut superlattices is plotted against t. The results for the spin-up channel, labeled as CoTt/Cut, and the spin-down channel, labeled as Colt/Cut, are shown for calculations with bulk disorder only, bulk plus spinindependent interface disorder, and bulk plus spin-dependent interface disorder. The lines are linear fits to the computed results. The interface resistances can be read from the intercepts. The bulk resistance is related to the slope. Both the interface resistance and the bulk resistance are affected by the types of disorder. It is found that both bulk disorder and spin-dependent interface disorder are necessary to explain experiment results [251. 1.3 Optimization at the Atomic Level Optimizations in physical systems are often performed by the two popular methods: simulated annealing [32] which is inspired by statistical mechanics, and genetic algorithms [33-34], which are inspired by evolution. Simulated annealing has been used in optimization problems such those in circuit design [32]. It is also used in problems such as structural optimization. For example, it has been used 15 0- Bulk disorder *-- -* Spin-indep interface disorder 3.0 o-- -o Spon-dep. interface disorder 2.0 1.0 Cot/CL 0.0 0 30 60 t (A) Figure 1.9: Current-perpendicular-to-plane resistivity of superlattices Co/Cu. The spin-up and spin-down channels are labeled CoTt/Cut and Colt/Cut, respectively. The resistivity are calculated for three cases: with only bulk disorder, bulk plus spin-independent interface disorder, and bulk plus spin-dependent interface disorder. The lines are linear fits to the computed results. The interface resistances can be read from the intercepts [25]. to find the most stable structure of a 7 Si atom cluster using the pseudo-potential density-functional theory [35]. The genetic algorithm has also been used in optimization of the molecular structure. For example, Deaven and Ho [36] performed a structural optimization by starting with a population of random molecules with 60 carbon atoms. Then the more stable, or "fit", molecules are cut into half and glue, or "mate", with each other to produce a new population of molecules. Genetic manipulations such as crossing over are performed in the process. As the simulation goes on, the most stable structure, the ideal icosahedral buckyball structure, is achieved. In Chapter 5, we turn to a new direction of research in which we optimize the atomic configuration of a nanostructure to obtain the maximum GMR. Performing optimization of physical properties in the atomic scale is an exciting new trend. 16 Franceschetti and Zunger performed the simulated annealing optimization to find the atomic configurations with the maximum bandgap in semiconducting alloys and superlattices [37]. To maximize the bandgap energy with simulated annealing (Figure 1.10), they started with a random atomic configuration and performed Monte Carlo moves in the configuration space, which is the space of all possible configurations. A move can be accepted with a probability even when the new configuration has a smaller bandgap. However, as the annealing goes on, the probability is reduced until the system finally rests in the maximum bandgap configuration. The maximum bandgap configuration for a A10.25Ga0.75As alloy is shown in Figure 1.11. It would be difficult to find this configuration using trialand-error. 190 b a 1.86 1.82 1.78 C M Al Gao 7,As alloy 1.74 0 2.000 4,000 6,000 8,000 10,000 12,000 Monte Carlo move Figure 1.10: Optimization of the bandgap energy of Al0.25Ga0.75As alloy by simulated annealing [37]. The bandgap energy fluctuates as the Monte Carlo routine progress. However, since there is a bias towards increasing bandgap and the amount of fluctuation is decreasing, the system finally arrive a maximum bandgap configuration. Although atomic level optimization has been performed, the optimization of physical properties other than the energy is new. Here one searches for an atomic 17 Figure 1.11: Maximum bandgap configuration of A10.25Ga0.75As [37]. This complicated configuration is unlikely to be found by trial-and-error. configuration that is, by definition, not the most stable configuration. However, with the non-equilibrium growth techniques, structures optimized for the desire physical property may become a reality. 1.4 Overview of the Rest of the Dissertation The rest of the dissertation goes as follows. In Chapter 2, we will develop a method to compute the conductivity in magnetic nanostructures. Both the band structure and the disorder will be included. In Chapter 3, we will discuss how to apply the model efficiently by rewriting the equations. In Chapter 4, we apply the method to calculate the CPP-GMR and compare the result with experiments. In Chapter 5, we will optimize the atomic configuration of a nanostructure to obtain the maximum GMR. In the final chapter, we will discuss the implications of this work and future directions. CHAPTER 2 MODELING TRANSPORT IN NANOSTRUCTURES In this chapter, we discuss how to model electron transport in magnetic nanostructures. The common approaches include semi-classical methods, simplified one or two band tight-binding model, multi-band tight-binding model, and densityfunctional calculations. Each model contains an increasing level of details and has its advantages and disadvantages. Semi-classical methods [16] can be used for studying large systems with complex geometry that other methods can not solve efficiently. The conductivity is often calculated with the Boltzmann equation. Bulk and interface scattering effects can be included as parameters in the Boltzmann equation, without the need for microscopic details. This approach valid when the length of smallest feature is much larger than the Fermi wavelengths, but it is often possible to extend it to smaller systems for qualitative study. However, when the feature length scale is only a few nanometers, the semi-classical method may become problematic. The simplified one or two band tight-binding model is often used to give qualitative understanding of the finite size effects and some of the band structure effects. Usually, the parameters in this model are obtained by reasonable estimations. The conductivity is usually obtained by the Kubo formula [38] or the Landauer formula [39-40]. This model is easy to solve, and sometimes it is even possible to obtain analytical solutions. In addition, the simplified band structure is helpful in understanding some effects that may be less obvious in more complicated 18 19 models. However, in real materials, such as iron or nickel, the band structure is too complicated for the results of this simple model to compare with experiments quantitatively. The multi-band tight-binding model [22-26] is used for quantitative comparison with experiments. Parameters are often fit from density-functional theory calculations and hence it captures the band structure. Details of atomic variation, impurity scattering effects [25-26] can be modeled accurately. Although this method is not as accurate as a full density-function theory calculation, it has the advantage of being much faster. It is also slower than the simple tight-binding model; however, as we will show in detail in Chapter 3, it is fast enough to perform large scale studies in inexpensive modern computer platforms such as a small Pentium III cluster. The density-functional theory is the most accurate method [27, 41-42], and unlike the above mentioned methods, less insight, or educated guess, is required in the calculations. The density-functional calculations can also provide information for the less demanding tight-binding calculations, and is useful for validating the results obtained from other studies. However, it demands much greater computer resources than tight-binding methods. Recently, large scale densityfunctional studies, such as the complex magnetic properties [43], are performed with the massively parallel processing computers. In the following sections, we first discuss the multi-band tight-binding Hamiltonian and the impurity averaged Green's function. Then we will explain the calculation of the conductivity using the Kubo formula and the vertex correction needed for current conservation. In the last section, we discuss how to obtain 20 tight-binding and impurity parameters corresponding to experimental situation and compare some of the theoretical results with experiments. 2.1 Modeling of the Band Structure In this section, we discuss the modeling of the band structure using a tightbinding Hamiltonian. Then we explain how to model impurities in samples and how to obtain the impurity averaged Green's function, which contains the band structure information. 2.1.1 Multi-band Tight-binding Hamiltonian The general form of a tight-binding Hamiltonian is given by H = H(r, r'; j, j'; a, o')cj cr'.',',, (2.1) r,rJJ 1,, where c . and crja are creation and annihilation operators of an orbital j with spin a at site r. Depending on the crystal structure and the details required in the model, the hoping terms, H(r, r'; j, j'; , 0'), are set to be zero when r and r' are third neighbors or further. However, in bcc crystals, the third neighbor is essential. In modeling the transition metals, it is often important to keep all the orbitals in the outermost shell: s, p, d. The coupling between the majority and minority spins H(r, r'; j, j'; T, 1) contributes to spin-flip scattering. For example, to include the spin-orbit coupling, a term gL - S is added to the Hamiltonian [44]. For transports in transition metals samples we consider, the spin-diffusion length is much larger than the layer thick- 21 ness [28-29], therefore spin-flip is neglected in most calculations. In this case, the two spin species are assumed to be independent and with decoupled Hamiltonians, H(T, T) and H(t, 1), and can be solved separately. Unless otherwise specified, the formulation is this work can be applied to the full Hamiltonian as well as the Hamiltonian for individual spins. 2.1.2 Modeling Impurities Scattering sources exist in all samples including in experiments of heterostructures. If we do not include these effects in our calculations, the resistivity computed would be at least one or two magnitudes less than those observed in experiments. There are three major source of scattering: disorder in the bulk of the sample, impurity atoms usually due to interface diffusion, and phonon scattering. Figure 2.1 shows a clean multilayer lattice and the same lattice with interface diffusion, where 10% of each monolayer diffuses into each of the two monolayers above and below. One can think of an impurity atom in a host as a new heterostructure. However, as we will explain in the next subsection, the advantage of treating it as an impurity is that simplification is often possible. For example, in an Co/Cu interface, we would treat the Co atoms diffused into the Cu layers as Co impurity in Cu host, and the Cu atoms diffused into the Co layers as Cu impurity in Co host. The impurity potential of an impurity B in an host A is defined as Au (r, r') = HB(r, r') - HA (r, r') (2.2) 22 where HB(r, r') and HA (r, r') are the hoping between r and r' for the perturbed and unperturbed system. The most important contribution to Au is usually the on-site term. As a perturbation, we found it sufficient to model this term as diagonal in all indices. In other words, Au(r, r'; j, j'; o, o') = Au(r, r; j, j; a, 0)6rr 6jj6'. (2.3) (a) (b) Figure 2.1: Illustration of a multilayer without and with interface diffusion. Without the interface diffusion (a), the multilayer has translational symmetry in the in-plane directions. In (b), 10% of each monolayer diffuse into each of the two monolayers above and below. By adding interface diffusion, the translational symmetry is broken. Two types of disorders are included in our model: the substitutional disorder which is often caused by interface diffusion and the structural disorder which are randomly located in the sample. To model the substitutional disorder in the 3d transition metals, we set the substitutional disorder parameter 23 AusUb(r, r; j, j; o, a) = 0 for j = s, px, py, pz because the most important difference between the tight-binding parameters in different elements is the on-site potential of the d-orbitals. The difference among the d-orbital parameters are small; therefore, it is often a good approximation to set AuSUb(r, r; j, j; o, a) = AUd(r, a) for j = dzy, dyr, dz, di, d2, such that two parameters, Aud(r, T) and AUd(r, 1), parameterize the impurity at site r. To model structural disorders, we use impurity potential of the form Au(r, r'; j, j'; o, o') = Austorrej'6,, (2.4) where Auq is a single parameter that characterizes the bulk of the sample. In this work, we do not model the phonon scattering. However, we can increase the structural disorder parameter for the bulk to get an idea of the effect of increasing temperature. 2.1.3 Impurity Averaged Green's Functions Once the Hamiltonian and the impurity parameters are given, we are ready to obtain the band structure information by either solving the eigenvalue problem or the Green's function. Here we explain the impurity average Green's function approach which is especially useful in treating systems with impurity scattering. The impurity average Green's function is the ensemble average of the Green's function for systems with different microscopic impurity configurations in the same macroscopic impurity distribution. In many systems, geometric symmetries are broken by the impurity atoms. For example, in the multilayer structure in Figure 24 2.1(b), the in-plane translational symmetries in the unperturbed system (Figure 2.1(a)) are broken by the impurities due to interlayer diffusion. However, the impurity probability distributions in each site within the same layer are the same. Therefore, after averaging over the impurity ensemble, the system recovers the translational symmetry. Also, the average of the Green's function is performed analytically. As a result, the impurity averaged Green's function approach can be used to study much larger systems than techniques that requires averaging over the impurity configuration numerically. We first consider the system unperturbed by impurities. The unperturbed, or bare, retarded/advanced Green's functions is given by ( - H ir)GO/A 1, (2.5) where the scalar w is the energy, rq is an infinitesimal small positive number, and 1 is the identity operator. The Green's function and the quantities obtained from the Green's function are w-dependent, in the following, the w label is suppressed for clarity. From the definition, the bare Green's functions satisfy (GR/A)t - GA/R. (2.6) Since the tight-binding parameters we use are real, the Hamiltonian is real' in the real space representation. In other words, HT = Ht = H, which leads to (G /^) = G/A (2.7) 'This is not true in the Fourier space 25 When impurity atoms are added to the system, one may try to compute the Green's function for system with impurities at their exact location. However, in most experiments, the exact location of impurity atoms are unknown, and even the probability distribution of the impurities are obtained by educated guesses. Moreover, the effect of the randomly distributed impurities are averaged out. Therefore, one has to compute the impurity averaged Green's function2, defined as (w - H - ER/A)GR/A = 1. (2.8) where ER/A is the self-energy. For the on-site, diagonal impurity potential described in the last section, the self-energy is given by ER/A (ri, r2) = v(ri)6rir2GR/A(ri, r2) (2.9) where the impurity parameter v(ri) = (Au2(ri)) (2.10) is the ensemble average of the square of the impurity potential. For convenience, we denote Equation 2.9 as ER/A = vGR/A. (2.11) 2For convenience, we use GR/A instead of (GR/A) as the symbol for the impurity averaged Green's function. 96 In this formula, we have neglected the cross-diagrams, which corresponds to fourth order terms in the impurity perturbation [45-46]. Since the impurity average Green's function and the self-energy depend on each other, we have to solve Equations 2.8 and 2.9 self-consistently. Other types of scattering effects, such as those by the bulk disorder, can be modeled within the same framework. We can also write the Green's functions in Dyson's equation as CR/A = Go/A + Go 0/ R/AGR/A. (2.12) The following properties follows from Equation 2.6: (GR/A)t = GAIR, (ER/A)t ZA/R. (2.13) In addition, in real space representation, (GR/A)T GR/A, (YjiRIA)T - ER/A. (2.14) 2.2 Conductivity Calculation In this section, we explain how to compute the conductivity using the Kubo formula which is the theory of linear response. We first review different current related operators in the tight-binding model, which will be then be used in the 27 Kubo formula. We will also discuss the vertex correction needed to give a current conserving conductivity. 2.2.1 Formal Aspects of Current Operators In this and the next subsections, we discuss discuss four operators related to the current: the total current, the divergence of current density, the current through a surface, and the current density operators. These operators will be useful for computing the conductivity and checking the current conservation. In this subsection, we will discuss these operators formally. In the next subsection, we express these operators in matrix form explicitly for an one dimensional Hamiltonian. Classically, the current density is defined as the charge density times the velocity. The total current, the divergence of current density, and the current through a surface can be obtained from integrating or differentiating the current density. In quantum mechanics, the current operator is defined using the Hamiltonian. However, in the tight-binding model, the current density is more difficult to define because of the discrete nature of the basis. The current density represented in a discrete basis, it is defined from the current from atomic sites r to r' divided by the cross-sectional area. In other words, there is only off diagonal element in the matrix representation. Unlike the classical current density, the current density operator not a localized in space. On the other hand, the total current, the divergence of current density, and the current operators are easily found from the charge density operators though the current conservation equations in Table 2.1. We will use the current operator to define the current name operator relation to Jr conservation equation total current J Er VrJr J =" = 1 Er Vrrpr divergence Vr Jr := JI(a+a.)-J"(x,) Vr - Jr = -APr current I Zr - Jr I E(r-ro).n Table 2.1: Summary of current density related operators in the tight-binding model. The names, symbols, relations to the current density operator Jr, and conservation equations used are shown in the table. In the tight-binding model, we define the current operators using the conservation equations, and use them to define the current density operator. In the table, Vr is the volume per atom, Pr is the charge density, P is the polarization operator, a, is the atomic spacing in the a-th direction. density operator, and use the operator for the total current and the divergence of the current density to check our definition. The relations between the operators discussed in this section are similar to those of the classical quantities; however, some details must be treated carefully. Table 2.1 summaries the current operators and how can they be obtained thru current conservation laws. The total current operator can be obtained [47] from the polarization operator, which is defined as P = e rcr. r (2.15) 29 The total current operator is given by OP att = i[H, P] = -ie Z(r - r')H(r, r')c crr. (2.16) r,r' Here, the Planck's constant divided by 27r, h, is set to equal to one. Since classically the total current is the volume integral of the current density, one may be tempted to define the current density operator Jr as proportional to Er,(1/Vr,)(r - r')[H(r, r')ctcr, - h.c.]. However, this definition violates current conservation because the divergence of this expression does not agree with the 0pr/Ot. We find the operator Vr Jr using the continuity equation Vr -Jr = &r (2.17) at In tight-binding model, the density is given by Pr = eccr /Vr, (2.18) where Vr is the volume occupied by the atom at site r. The Schr6dinger equation implies Vr Jr = -i[H,,pr] -ie = [HctCr - CtCr H] (2.19) Vr r r 30 -- [H(r', r)c,cr - rc H(r, r')]. (2.20) r r' It is useful to note that in matrix representation the only nonzero element of ctCr' is 1 at (r, r'). In other words, (Ctcr')(rl, r2) = 6rrl 6r'r2. Therefore, we can write the element of the divergence of the current explicitly as (Vr Jr) (ri, r2) = (-ie/Vr)[H(ri, r)6rr2 - rir H(r, r2)] (2.21) The operator for a current through a plane containing point ro and normal to a unit vector n can be found by a I -= - E PrVr (r-ro)-n - -ie E H(ri, r2)rCr2 , Z Ccr] rir2 (r-ro)-n (r-ro).n =--ie ( E + E ) [H(ri, r)cc - h.c.] \(r-ro).n (r-ro).n We have used the fact that the second term in the fourth line is zero due to the symmetry between r and r1. Physically, the second term is zero because it contains only the hoping term that does not leave one side of the plane. The plane current operator obtained at the end contains only the hopping terms which bring current from one side of the plane to another. 31 Unlike the operators shown in the first column of Table 2.1, there is no simple equation that leads to the current density operator. In the continuous limit, the current density at r is given by the current though a infinitesimal surface divided by the area. However, in the tight-binding model, we have to discretize the length scale to the atomic scale. As an example, in one dimension the current density operator in second quantized form is given by J. = -ia/Vt E [H(x1,x2)c cx2 - h.c.]. (2.23) X1 X>X2 One can check that the divergence, defined in one dimension as (Jx+1 - Jx)/a, agrees with Equation 2.27, and Ex JxV equals the total current J. This definition can be extended to higher dimensional lattices. For example, in a rectangular lattice with only nearest neighbors, the two component of the current density is given by J,1 = -ia1/V [H(ri, r2)cr,2 - h.c.], X1 X>X2,y1 =Y2=Y -ia2/Vr [H(ri, r2) Icr2 - h.c.], (2.24) Y1 Y>2,X1=X2=X where a, and a2 are the lattice spacings in the x- and y-directions respectively. 2.2.2 Current Operators in One Dimension In this subsection, we consider an one dimensional Hamiltonian and express the corresponding current operators in matrix form explicitly. We can better un- 32 derstand the properties of the current operators by studying this one dimensional example. Consider a one dimensional Hamiltonian, H, with up to second neighbor hoping terms3. In matrix form, the Hamiltonian can be written as The Hamiltonian in matrix from is given by ( H(5, 3) H(5, 4) H(5, 5) H(5, 6) H(5, 7) H(6, 4) H(6, 5) H(6, 6) H(6, 7) H(6, 8) H(7, 5) H(7, 6) H(7. 7) H(7, 8) H(71 9) (2.25) From Equation 2.16, the total current operator J can be found by -iea times the matrix 2H(5, 3) H(5, 4) 2H(6,4) 0 H(6,5) 2H(7,5) -H(5, 6) 0 H(7,6) -2H(5, 7) -H(6,7) 0 -2H(6,8) -H(7,8) -2H(7,9) where a is the lattice constant. 3The model with only first nearest neighbor is too simple to illustrate some of the properties. / 33 The matrix form of the current density operator at, for example, x = 6, can be obtained from Equation 2.23 as -H(4, 6) = -iea -H(5, 6) -H(5, 7) V H(6, 4) H(6, 5) H(7, 5) (2.26) where a is the lattice constant, and V is the volume occupied by a site. using this definition, one can see that the total current operator satisfies J = Ex JxV. We also note that in one dimension, the current operator is equals the current density times the area V/a. In other words, the current operator at x = 6 contain only the hoping terms which bring the current from one side of the plane x = 6to the other side. By Equation 2.26, the divergence of the current density operator at x = 6 is given by d - Jx dx J x=6 J7 A a 34 I H(4,6) H(5,6) -H(6, 4) -H(6, 5) 0 -H(6, 7) -H(6,8) I - H(7,6) H(7,7) / (2.27) This equation can be rewritten in the following form: H(4,6) H(5,6) 0 0 H(6, 6) H(7,6) H(8,6) ( d -Jx dx x=6 V 0 0 'I = e 35 0 0 +V H(6, 4) H(6, 5) H(6, 6) H(6, 7) H(6, 8) 0 0 (2.28) which the same as Equation 2.20 in one dimension. If we sandwich this expression by two Green's functions or wavevectors, the first and the second terms of the expression can operate on the left and right respectively like the Hamiltonian. For example, the second matrix operated on the retarded Green's function on the right can be evaluated using Equation 2.8: E( - H(x, x1) - ER(x, xi))GR(Xi, XI) = 6(x, X'). (2.29) X1 This property of divergence of the current density is useful later in the discussion of current conservation. 36 2.2.3 Kubo Formula In the linear response theory, the non-local conductivity tensor is defined as the response function in the equation Jr = EZ-rrErVr,, (2.30) r/ where Jr is the current density response at r, Er' is the applied electric field at r', and Vr is the volume per atom. The impurity averaged conductivity tensor can be computed by the Kubo formula : Urr, = 1 Tr [2Jr,G AFARGR - Jr,G RFRRG R- jrGA A G . 47rrrr r (2.31) In these equations, the current density plus vertex correction JrAR is given by the Dyson's equations r Jr + Kr (2.32) = Jr + vGA rAR GR (2.33) =Jr + VGA JGR +vGA vGA JrGR G +.... (2.34) The other vertex corrections KA, K R, and KA^ are obtained by replacing the corresponding indices in the equations above. The Kubo formula is represented in the Feynman diagram is Figure 2.2. Since, Kr is diagonal and 1r is not, in practice, it is more convenient to iterate on KrR instead of FR: KR vGAFARG (2.35) = vGA JGR +VGKRGR. (2.36) We also define 0= >r Orr' Vr' r/ 1 Tr [2JG^ ARGR - JGR r RRGR - JGA rAGA . (2.37) For a constant electric field, Jr = o-rE. Therefore -r is often compared to experiments. It can be shown that this conductivity is still correct even in non-uniform electric field because the result depends only on the potential drop across the sample. As seen in the bubble diagrams in Figure 2.2, the conductivity can be found by first dressing up the vertex in either current operator and then multiplying it by the other current operator, followed by taking the trace. Analogous to J, we define f^-Er Vr and KRA r rRAgr, which gives r RA = J KRA (2.38) = J+vG RT RAG (2.39) = J +vGRJG+A vGR vGRJGA G+A .... (2.40) 38 Similarly, the iteration equation for KRA is KRA = vGR^RAG = vG+JG vGK GA. (2.41) (2.42) Therefore, we also have Or= 1 Tr [2r RAGAJrGR - rRRG RJrGR - rAAGAJrGA 47r (2.43) As we will discuss in the next chapter, the two sets of equations, Equations 2.322.37 and 2.38-2.43, have different computational costs in different situations. In Chapter 3, we will explain which set of equations should be used in currentperpendicular-to-plane and current-in-plane calculations. (a) (b) A R A R -x-- A R + +X~ + -X- Figure 2.2: Feynman diagrams for the average conductivity. (a) The open bubble, which corresponds to the conductivity without the vertex correction. (b) The conductivity contribution from the vertex correction. 39 2.2.4 Current Conservation For a system without impurity average, the conductivity is given by the the Kubo formula without the vertex correction, or the open bubble (Figure 2.2(a)). However, when averaged over the impurity configurations, the open bubble does not satisfies current conservation. This is because using the impurity averaged Green's function to compute the conductivity in the open bubble, one neglects the correlation terms such as those in Figure 2.2(b). These terms are required to give a conserving current. In the following, we show the terms in the vertex correction shown in Equation 2.34 will give a conserving current. In other words, we will show that the conductivity obtained from Equation 2.37 satisfies Vr - Jr = Vr - UrE = 0. (2.44) Since the electric field is arbitrary, we need to show Vr - 0r = 0. (2.45) Since the current Jr, does not depend on r, we can first find the divergence of the vertex GAFrGR, and then multiply the divergence with the Jr,, followed by taking the trace. The divergence of the zeroth order term of the AR vertex in Equation 2.34 is Vr - (GA JrGR) = GA (Vr - Jr) GR. (2.46) 40 In the one dimensional case, this is simply (G Jx+1GR - GJxGR) = A(Jx +1 - X )GR. a By current conservation (Equation 2.19), we have GA(Vr - J,)GR = (-ieVr)[(GAH) ccrGR -Gccr (HGR)] = (-ie/Vr)[(-1 + GAW - GA EA)ctCrG R - GAc Cr(-1 +wG R = (-ieVr)[G rcr - CtrCrG - GA ArCrGR + c G RCr G (2.47) ERG R) (2.48) Here, we have used the Schr6dinger equation (2.8) to evaluate HGR/A as (2.49) For the divergence of the first order term of the AR vertex (vertex with one rung) in Equation 2.34, we have Vr - (G AvG JrGR ) GR) = GA v(GA(Vr Jr)GR GR = (-ie/Vr ) GA vGCtCr G R - GA VCtCrG R G R -GA vGAE crG GR +G A vGCCrR GR G R (2.50) HGR/A = -1 + wGR/A - ER/A GR/A 41 The first term of the above equation can be rewritten using the following equation: (vGA ccr)(ri,r2) v(ri)6rir2(G'ctcr)(ri,r2) = v(ri)6r1,r2GA (ri, r)6r,r2 =EA (ri, r)6r,r2 (EActrr)(ri,r2). (2.51) We can see the first term of Equation 2.50 cancels the third term of Equation 2.48. Similarly, the second term of Equation 2.50 cancels the last term of Equation 2.48. The same argument can be applied for higher order terms to show the last two terms of the divergence of the (n + 1)-th order vertex cancel the last two terms of the divergence of the n-th order vertex. As a result, the only terms left of the divergence of the vertex are the first two terms of Equation 2.48. The divergence of the conductivity in Equation 2.37 is proportional to the sum of the following terms 2Tr(JGAVr - 4A)G) (-2ie/Vr )[i - r21] -Tr(JGRV, (FRR)GR) = (ie-Vr)[, KR -Tr(JGAVrA (AA)GA) = (ie/Vr)[t - ], (2.52) where ^ - Tr(PJcrGR/^) = Tr(crcrG R/J), RIA _ Tr(JGR/^Accr), 42 (2.53) In the real space representation, HT = H, which implies (ccrG RAj)T - JG R/ACt Cr . (2.54) Since transpose of a matrix does not affect its trace, we have R/A _ R/A (2.55) Unlike the continuous case [46], the divergences of the individual dressed bubbles in the tight-binding model are non-zero. However, they cancel each other so that the divergence of Ur obtained from the Kubo formula is zero. Therefore the vertex corrected conductivity satisfies current conservation. In this work, we only use on-site scattering potentials, the self-energy ER/A is also on-site. Therefore, EA/Rctc_ = 4Cr .rt . (2.56) This means that in Equation 2.48 the last two term for divergence of the AA and RR vertex cancels each other out, and no vertex correction is needed for the AA and RR vertex. The Kubo formula for on-site scattering potentials can be simplified as 1 -r = Tr [2JGA(Jr + KAR^)GR - JGRJrGR - JGAJrG^ , (2.57) 43 When the conductivity is calculated in nanostructures, for example in the current-perpendicular-to-plane (CPP) magnetic multi-layers, we find the contribution from the vertex corrections can be as large as 20% of the total conductivity. Therefore, the vertex correction should not be neglected. 2.3 Realistic Parameters In the earlier part of this chapter, we explained how to obtain the conductivity once the Hamiltonian and the impurity potential are given. In this section, we explain how to obtain the tight-binding parameters and the impurity potential, which are essential to modeling the electron transport accurately. We will also compare some theoretical results with experiments. 2.3.1 Tight-Binding Parameters The tight-binding parameters of an element is obtained from fits to the bulk band structure of density-functional calculations [48]. To do this, we compute the eigenvalues for a set of k-points in the irreducible zone with the density-functional theory, and also with the tight-binding model. The tight-binding parameter are then adjusted to minimize the error ZEZ LTB ( WPFT (k)2, (2.58) kj where ,TB(k), and wPFT(k) are respectively the tight-binding theory eigenvalue and density-functional theory eigenvalue of the j-th orbital of a k-vector. As shown 44 by Papaconstopoulos [48], tight-binding parameters obtained in this way give the correct band structure for bulk systems. For heterostructures, we do not obtain the tight-binding parameters by fitting to the band structure of the heterostructure because the density-functional theory for complicated systems can be expensive computationally. Since our goal is to model the system efficiently but still give reliable predictions, we use a less costly approach. The band structures for 3d transition metals are very similar, except for the offset of the d-bands, which are different in each metal. To model 3d transition metal heterostructures, we use the bulk parameters in the bulk and the interface. The on-site energy of each metal is shifted so that the Fermi-level level is zero. In this case, the Fermi level for different metal aligns. The hoping parameters between the atoms of two difference element are set as the average of their bulk values. 2.3.2 Impurity Potential for Substitutional Disorder As explained in Section 2.1.2, the substitutional disorder is caused by impurity atoms usually from interface diffusion. The impurity potential for the atom can be obtained from density-functional fits. To the first order, one may use the difference in the bulk on-site parameters between the impurity and the host as the substitutional disorder potential, ansub. However, since in transition metals the on-site energy of the impurity B are often affected the host A, it is more accurate to find the impurity potential by fitting the tight-binding model to the density functional theory for a small unit cell such as A3B1. Again, the error between the tight-binding and density-functional eigenvalues given in Equation 45 2.58 is minimized. Since in the transition metals the most important difference between different elements is the d-orbital on-site energy, all the parameters except for the d-orbital on-site energy of the impurity can be set as the parameter of the host obtained from the bulk fit. As a result, the fitting for A3B1 is just an oneparameter fit for each spin, where the only varying parameter is the d-orbital on-site energy. The validity of the parameters are then confirmed in small unit cells by comparing the density of states computed by the tight-binding theory and the densityfunctional theory. For example, parameters for Fe as impurities in fee Ni host are obtained from fits to fec Ni3Fe band structures. As a check, parameters obtained from the fit are used to calculate the tight-binding band structure of fee Ni3lFe, which agrees with the band structure obtained from the density-functional theory. This indicates that the model works well at least when impurity concentration is less than 25%. When the impurity concentration is large, or when the difference in hoping between the host and the impurity is large, the error is expected to increase. The other parameter in interface diffusion is the density of impurity atoms in the few atomic layers next to the impurity. Unfortunately this parameter depends the material and methods used to grow the sample in experiments. Only in a few transport experiments the sample is characterized well enough to give the details of the interfaces. Even in those case when the interfaces are well characterized, such as those by Cyrille et al. [6-7], only the long length scale (~10 nm) interface roughness is measured. Within our model, the interface roughness at this length scale should be considered in the geometry because the Fermi wavelength is an 46 order of magnitude smaller. On the other hand, the atomic scale interfaces roughness is related to the density of impurity atoms. In Chapter 4, we will explain how to set the density of the impurity atoms. 2.3.3 Impurity Potential for Structural Disorder The impurity potential for the bulk structural disorder, Aust, can be obtained by comparing the calculation of the resistivity and the resistivity for bulk samples. Usually, the samples in heterostructures are not as clean as a bulk sample. Therefore, we use the room temperature resistivity in bulk samples as a guide to setting the bulk disorder parameter. 2.3.4 Extension to Magnetic Alloys The method described in this section can also be extended to model transition magnet alloys, where the magnetic moment and hence the band structure vary with the composition. In these materials, the localized d electrons are responsible for the magnetism. We can assume that the most important changes due to alloying are contained in the on-site energies, H(r, r, d, d, T, T) and H(r, r, d, d, 1, 1), of the spin-up (majority) and spin-down (minority) d electrons. In other words, the alloy is assumed to have site-independent hoping parameters and s and p on-site energies the same as the host metal, and site-dependent d on-site energies. This assumption will be justified by the agreement with local density calculations of supercells in different impurity concentrations. The major contributions to the on-site energies come from the atomic core potential plus the Coulomb energy due to the opposite spin. In the itinerant 4- electron model, it is more convenient to work with the total number of spin-up and spin-down d electrons, N and Nd, at each site r. Therefore, we rewrite the parameters in the form of H(r,r,d,d, ,I) = U+UNd H(r, r, d, d, ,, I) = UO + UxN1, (2.59) where Ur? and Ux are the on-site parameter and effective Coulomb energy per pair at site r. Both N1 and N' in Equation 2.59 are obtained from the band structure, which in turn depends on H(r, r, d, d, 4, 1) and H(r, r, d, d, 1, t). Thus, Equation 2.59 has to be solved self-consistently with the alloy band structure calculated using Hamiltonian of Equation 2.1. The alloy band structure is computed with the coherent potential approximation (CPA) [49-50] is used, which is summarized in Appendix B. All parameters for the host atoms can be found from fits to the local density calculations of the pure metal [48]. The only two parameters left, UmP and Uxn of impurity atoms, are obtained from fits to the local density calculations of supercells composed the two types of atoms. In Figure 2.3, we plotted the Slater-Pauling curve, which shows the magnetic moment against the average number of electrons, for 3d magnetic alloys. Both experimental data and tight-binding CPA calculations are shown. As we can see, the theory agrees not only with the main curve, but also the NiCr branch as well. The NiCr curve is very sensitive to the details of the band structure. In this simple model, only two parameters in addition to the parameter of the host are 48 included; therefore, there is some deviation in the NiCr curve. The agreement indicates that the tight-binding model described here is suitable for obtaining the magnetic properties and band structure properties of transition magnetic metal heterostructures. E 0 I CL a> E 0 E CM C E 6 7 8 9 10 11 number of valence electrons per atom Figure 2.3: The magnetic moment of alloys per atom versus the number of valence electrons per atom. The lines are the experimental data of FeV, NiFe, NiCu, and NiCr alloys as labeled. The symbols are from the CPA calculations of FeCr (triangles), NiCr (circles), NiCu (squares), and NiCr (crosses). The CPA calculation agrees well with the experiment in the main branch (NiFe and NiCu), and agrees reasonably well with the NiCr branch. The dip in the experiment curve in NiFe corresponds to the bcc-fcc phase transition. FeV NiFe NiCu A I I sa NiX 21 1 CHAPTER 3 COMPUTATION PROCEDURE In the last chapter, we discussed the formal aspects of modeling magnetic nanostructures and calculating the conductivity. We also discussed how to obtain parameters so that the models give realistic results that can be used to compare with experiments. However, if one writes computer codes directly from these equations, the resulting codes will not be fast enough for large scale studies. In particular, many of the equations involve O(N3) operations such as the matrix multiplication of NM x NM dense matrix, where N is the number of atoms in a unit cell, and M is the number of orbitals per site. In order to speed up the calculation to tackle the tasks such as optimization and studying large systems, we have to reformulate the problem with special consideration about the computation speed. In this chapter, we discuss the computation procedure of the conductivity calculation. In the first section, we discuss the current-perpendicular-to-plane (CPP) conductivity computation, which demonstrates how to find the impurity averaged Green's function, solve for the vertex correction, and compute the conductivity. Although we use the CPP geometry as an example, the application to general geometry are similar. Most importantly, the tools that are developed here can be used in more general cases and the size-dependence of the computational cost are the same. As we will show, the size-dependence for computing the self-energy, the vertex correction, and the conductivity are O(N)., O(N2), and O(N), respectively. 49 50 This is an O(N) saving compared to total computation costs in the last Chapter. As a result, we can study more systems with a larger size. As we will see in Chapter 5, the speed of the code is crucial to optimization. In the second section, we will discuss the current-in-plane (CIP) geometry, which is also commonly studied in experiments. We will show that since the Hamiltonian has translational symmetry in the in-plane direction, the vertex correction is identical to zero. Therefore the CIP conductivity is easier to computer. However, the size-dependence of the overall computation remains O(N2). In the third section, we briefly discuss the calculation of the conductivity tensor and the current-at-an-angle-to-plane (CAP) geometry, which will only require tools developed in the first two sections. In the last section of the chapter, we explain how we implement the algorithms discussed in this chapter. We shall use some new C++ matrix libraries, which will call the FORTRAN libraries BLAS (Basic Linear Algebra Subprograms) [51] and LAPACK [52] for the actual computation. We will also use an in-house blockmatrix library to implement most algorithms, and explain the advantage of using this new tools for coding. 3.1 Current-Perpendicular-to-Plane Conductivity We discuss the computational aspects of the current-perpendicular-to-plane (CPP) conductivity calculation in this section. Let us consider a periodic multilayer system that has N mono-atomic layers in each period and can be generated by a unit cell of N atoms. We label the growth direction as z and the in-plane directions as x and y. Since in the unit cell each atom belongs to a different mono- .31 layer, we label the an atom in the unit cell as z = 1, 2, . N. An example of this is a multilayer grown in the fcc(111) direction. Each elements of a matrix is an m x m block matrix, where m is the number of orbitals per site. The current is in the z-direction. The equations in this section are intended for direct implementation, therefore, the computational cost of major equations are given. 3.1.1 Impurity Averaged Green's Function Using the model discussed in Chapter 2, the self-energy is given by ER(ri, r2) = v(ri)GR(ri, ri)6(ri, r2) (3.1) We only need to calculated the real space function for sites within the unit cell. Since E is diagonal, it can be stored in a vector and the computation cost is N. In real space, the Green's function is given by GR = (w - H - ER (3.2) However, since a periodic system is also infinite, we have to calculate the Green's function for the unit cell in the Fourier space: (3.3) GQ =Ro Hk - ER)-- 52 where k is a wave-vector in the first Brillouin zone, Hk and ER are the Hamiltonian and self-energy in the k-space. Since the self-energy is on-site, we have ER = ER in the unit cell. The Hamiltonian in k-space can be found by Hk = ZR H(O, R)exp(-ik - R), where R is the position vector of neighbor unit cells. In the following, we consider up to second neighbor hoping in real space. In the fec(111) lattice, up to the second nearest neighbors are included in the three monolayers described by z = 0, 1. Therefore, when we work in the k-space representation, the Hamiltonian is tridiagonal with periodic boundary condition. In Appendix A, we explain how the inverse can be found efficiently. Since Hk is tridiagonal with periodic boundary condition and E is diagonal, the cost for inverse is O(N2). However, since only diagonal elements of G are needed to find the self-energy, the cost can be cut to O(N) in this step. The real space Green's function within the unit cell is then given by 1 R GR = E G (3.4) Nk k The advanced self-energy and Green's function can be found by taking the hermitian conjugates of the results. The Equations 3.1 and 3.4 are solved self-consistently by iteration. If the system is finite instead, then we can work in the real space and the inverse algorithm for fixed boundary condition in Appendix A can be used; however, one has to consider the lead with special care. 53 3.1.2 Vertex Correction To find the conductivity, we have to compute not only the diagonal of the Green's function, but also the full Green's function. In Appendix A, we explain how to compute the Green's function in O(N2) for both periodic and fixed boundary conditions. In the k-space, the current density operator at site z in the unit cell is given by' Jzk(zl, z2) = (-ia3/V)Hk(zl, z2) if z1 = z, z2 - z + 1 = (ia3/V)Hk(z, z2) if z1 = z + 1, z2 = z = 0 otherwise, (3.5) where a3 is the monolayer distance. The total current is Jk JzkV, (3.6) which has 2N elements. As discussed in Subsection 2.2.3, both sets of equations, Equations 2.32-2.37 and 2.38-2.43, can be used to compute the conductivity. When we compute the vertex correction, both Equations 2.32 and 2.38 have the same computational costs. To see this, we notice from the definition of KAR that it is non-local. Therefore, both KzR and KRA have the same number of elements and are computed in the 'Here the coordinates are cyclic, for example, z = N + 1 refers to the same site as z = 1. 54 same number of steps. Since J, is local, the cost for computing Tr[GFGJz] in Equation 2.43 is O(N), while the cost of Tr[GJGrz] in Equation 2.37 is O(N2). The vertex correction KRA is given by the Fourier transform of the Dyson's equation (Equation 2.36): KRA(z1, zi) = > v(zi)G'(zi, z4)JA(z4, z5)GA(z5, z1) k Z4,Z5 + 5 v(zi)Gk(zi, z4)KRA (z4, z4)Gk(z4, zi) Z4 v(zi) Ey[G (zi, z4)JA(z4, z4 + 1)GA(z4 + 1, z1) k Z4 +Gk(zI, Z4 + 1)JA(z4 + 1, z)Gk(z4, zi)] +v(zi) E[G (zi, z4)KRA (z4, z4)Gj (z4, zi)] Z4 (3.7) Again, Kk"A = KRA in the unit cell because KRA is diagonal. The first and second term cost O(N2), however, they are independent of the iteration. In the last term, for each zj, we need to do a sum over z4, thus, each iteration in this loop costs O(N2). Since this step is inside the iteration loop, it is the most time-consuming part of the whole calculation. The cost for computing GR to start the loop is also O(N2). 3.1.3 Conductivity The conductivity is given by Fourier transforming Equation 2.43 as -Z = I E Tr [2JzkG (Jk + KRA)GA - JAGR JkGR - JkGAJkG A 4rkkk k - zk k k .(3.8) ,55 As explained in Chapter 2, o, is independent of z, and hence we can choose any z the evaluate the conductivity. One should note that it only take O(N) to evaluate the traces. For example, the first two terms can be written as 'Tr [GkJzkG J = (G JzkGk)(z' + 1,z')Jk(z',z' + 1) ZI +(GA JAkGR)(z', z' + 1)Jk(z' + 1, z'), r [G A JAiG KkR^ = (G A JAGR) (Z', Z') KkA (z', Z'), and (3.9) (3.10) where (Gd JeG2/A)(z' + 1, z') (Gk JzkG2/^)(z', z' + 1) = (Gk JzkGR/^)(z', z') = G Z + 1, z)Jzk(z, z + 1)G RA(z + 1, z') +G + 1, z + 1)Jzk(z + 1, z)Gk/A(z, z'), Gj(z', z)JA (z, z + 1)G /A(z + 1, z' + 1) +G (Z', z + 1)Jzk(z + 1, z)GR/A (z, z' + 1), Gjd(z', z)J(z, z + 1)G2/^(z + 1, z') +G (z', z + 1)Jzk(z + 1, z)G /A (z, z'). (3.11) Note that there is only one sum in each trace, therefore, the cost for each trace is O(N). Further simplifications are possible by using the symmetry of the matrices, 56 for example, Tr[(G^JGR')(z, z'+ l)J(z' + 1,j)] = Tr{[J(z', z' + l)(GAJGR )(ZI + 1, z')]t} = (Tr{(GAJGR)(zI + 1, z')J(z', z' + 1)})*. (3.12) Therefore, Tr[GAJzGRJ] = 2 E Rj{(GAJzGR)(Z' + 1, z')J(z', Z' + 1)}. (3.13) 3.1.4 Diagonal Impurity Potential Most of the discussions in this work allow the impurity potential to be offdiagonal in the orbital indices. However, for the diagonal impurity potential, v(z)(m, n) = v(z)(m)6(m, n), further simplification that reduces the cost of the second loop to negligible is are possible. In this case, since KRA is also diagonal in band index, we denote KRA(z)(m) = KRA(z, z)(m, M). In the following, we will show that the last term in Equation 3.7 can be written in the form of a matrix-vector multiplication. Hence we can group the KR^ terms on both sides of Equation 3.7 and write this Dyson's equation in the form of (1 - q)K RA = , (3-14 where q is a matrix and p is a vector. Equation 3.14 is a set of linear equations and can be solved without iterations. To derive this equation, we notice that for each zi, KRA(zi)(m) = p(zi)(m) + v(z(m) G'(zi, z2)(m, n)K(z2)(n)GA(z2, zi)(n, m) Nk k Z2kk = p(zi)(m)+ q(zi, z2)(m, n)KRA (z2)(n), (3.15) Z2 where pAzi)(m) E v(zi)(m)E(G JkG')(zizi)(m,m) p~i( )Nk k k k _v(zi)(m) (lA = E( [Gk(zi, z2)Jk(Z2, z2 + I)Gk(Z2 + 1, zi) Nk k Z2 +G R(z, z2 + 1)Jk(z2 + 1, z2)GA(z2, zi)](m, m) v v(zi) (mr) (lA = >3>3[Gk(ziz2)Jk(z2,Z2+1)G (z2+1,zl)](m,m)+h.c. INk k Z2 _2v (zi) (m) ZIA ER{[Gk(ziz2)Jk(z2, 2 + 1)Gk(z2 + 1, zl)](m, m)}. Nk k Z2 (3.16) and q(z,z2)(m,n) v(z m)Z2)(m, n)Gk(Z2, zi)(n, m) = v(z m)(m, n) (3.17) 56 Since the second term of p is the conjugate of the first term, the only term in p we need to compute is the first term. It can be calculated as follows [G'(zi, z2)Jk(z2, z2 + 1)G(z2 + 1, zi)](m, m) = [G'(zi, z2) Jk(z2, z2 + 1)](m, m')G' (z2 + 1, zi)(m', M), (3.18) which does not involve any sum over positions. Since q is symmetric, we only need to calculate half of the element, and we can use the symmetric matrix solver for this problem. These two cost saving steps are important because they are O(N2) calculations. When impurity potential is not diagonal in the orbital indices, KR^ and G do not commute and it is very costly to solve the linear equations, therefore KRA is solved by iteration. However, for the impurity potential used here, it is much faster to Equation 3.14 directly. which is a set of real and symmetric linear equations and can be solved with a cost of N2M2/2. The cost for computing p and q are proportional to NkN2M3 and NkN2M2, where Nk is the number of k-points. Therefore, most time is spent in evaluating p. 3.1.5 Total Computational Cost The computational cost is important to the large scale we discuss in the following Chapters, therefore it is given here. If N, are the number of iterations in the self-energy loop respectively, and Nk is the number of k-points, then the total cost for a conductivity calculation is roughly (4M3)Nk(6N2 + (19N, + 27)N) floating 59 point operations (flop), where the factor 4M3 is the cost for a complex 9 x 9 block matrix multiplication. As an estimate, we take Nk ~ 1700, M = 9, N = 50, N, = 10, the total cost is on the order of 10" flop. On a machine capable of about 10' flop per second (1 GFlops), this should take a few minutes. In practice, we performed the calculations a computer cluster with ten 866MHz Pentium III CPUs, it takes about 5 to 10 minutes for a calculation when N = 50. The algorithms presented in this chapter are CPU intensive. The memory requirement is O(N2). For the computers used in this work, unless N is larger than about 500, one does not need to rearrange terms to optimize the memory usage. 3.2 Current-In-Plane Conductivity The current-in-plane (CIP) conductivity is a special case where the Hamiltonian has a translational symmetry in the direction of the current. Let us consider the same fcc(111) lattice discussed in the last section and denote the direction of the current as x. The problem is greatly simplified because the vertex correction KRA is zero. This result can be shown by studying the first correction term of Equation 2.40: (vGR JGA)(x, x) = v(x)[E GR(X, X + j)J(x + j, x j + 1)GA(x + j + 1, X) + GR(X, X +j + )J(X+j + 1, x + j)G A(X + jX). (3.19) 60 In the real space representation, the Hamiltonian is real and has translational symmetry, hence GR and G^ are symmetric, and J is antisymmetric. The second sum above can be rewritten as - GR(x + j + 1, r)J(x + j, x + j + 1)GA(X, X + j) EGR(x, x - j - 1)J(x -j - 1, x - j)GA(x - j, x), (3.20) where in the last line we have used the translational symmetry separately for GR, J, and GA. By substitute j' = j + 1, we can see the two sums in Equation 3.19 cancel exactly. As a result, all the correction terms to fRA in Equation 2.40 are also zero. Since the vertex correction is zero, we can compute the conductivity corresponding to the open bubble by 1 o2 = 41 ZTr[2JkGA JXkGR - JkGR JAkGR - JkGA JAkGA]. (3.21) Since J~k is a tridiagonal matrix in the unit cell, it has 3N elements. Thus, the computational cost for each trace is O(N2). 3.3 Conductivity Tensor As explained in Chapter 2, the current density operator is not obvious to construct. It is especially difficult when the direction of current is at an angle to the crystal axes. Therefore, in computing the current for systems such as the current- 61 at-an-angle-to-plane (CAP) system, it is easier to compute the conductivity tensor and use the tensor equation J = uE. In a cubic bulk system, o,1 = 6,0o". However, in nanostructures, the symmetry is usually broken, therefore a'0 are in general different. The techniques for computing current with vertex correction is useful for computing the conductivity tensor. 3.4 Coding As shown in the above, it is natural to formulate our problem in block matrices and hence natural to write our code using block matrix object classes. We have to implement the block matrix algorithms such as inverse and products discussed in previous sections to take advantage of the sparsity of the matrices. Therefore, it is not enough to just have a block matrix class and common algorithms. The block matrix class should help us develop new block matrix algorithm efficiently. Our goal is to translate mathematical equations involving block matrices into programming lines directly, while having performance close to Fortran codes. In this section, we discuss how this is achieved with C++. First we use the Template Numerical Toolkit package to handle the operations of all dense matrices. To increase the speed, we link critical subroutines, for example, the matrix multiplication or inverse operation for dense blocks, to the Fortran/Assembler coded native binaries BLAS [51] and LAPACK [52]. Next we develop a block matrix class that can use the TNT matrix class as elements. Using the new block matrix class, we develop codes to implement the algorithms described by previous sections. 62 In addition to these steps, our code is parallelized with PVM (Parallel Virtual Machine) at the level of k-points to use the all available computing resources of a heterogeneous network. The resulting code should scale well with the number of processors. 3.4.1 Template Numerical Toolkit (TNT) The template [53] feature in C++ language allows writing a single code for object classes or algorithms for different storage types. For example, the subroutine template T f(T x) { return x*x; } defines a function that can be called by f (x) whether the data type of x is int or complex< double >. The compiler generates suitable binaries to match for the data type. When we started developing the first version our code, the ANSI C++ standard was not formally established yet. However, there were two attractive experimental C++ packages that handle matrices well by making use of the template features: the Template Numerical Toolkit (TNT) by R. Pozo [54-55] and Blitz++ by T. Veldhuizen [56-57]. Although Blitz++ has great potential, it used advanced features such as expression templates, which were not supported by the compiler on 63 our platforms at the time we started coding.' On the other hand, the TNT package could be compiled and linked with Fortran programs on our major platforms: Solaris and Linux. Consequently, we chose the Template Numerical Toolkit. The TNT object we use most are the FortranMatrix object, which stores a matrix using Fortran (column) convention, and the Vector object. They allow Fortran (one-based) indices such as A (2,3). For example, the following lines define a 2 x 2 matrix A, vectors x and b. FortranMatrix< int > A(2,2, " 1 2 1 4 "); Vector< int > x(2, " 1 2 "), b(2); b(1)=3; b(2)=4; The package also takes advantage of the function/operator overload feature of C++ to make the matrix and vector operations intuitive. Basic operations like =, +, *, and << are overloaded. Continued from the above example, the next two lines computes and prints out the vector c = 3A2X + b. Vector< int > c=3*A*A*x+b; cout << c << endl; Notice that the compiler links the three multiplications above to scalar-matrix multiplication, matrix-matrix multiplication, and matrix-vector multiplication subroutines in order. Without operator overloading, the line that computes c would look like Vector< int > c(vecadd(vecmult(matmult(scalmult(3,A) ,A) ,x) ,b)); 2According to the Blitz++ website [571, Blitz++ can now be compiled by compilers like gcc and KAI C++. 64 3.4.2 Fortran performance with BLAS and LAPACK To achieve the performance close to that of a Fortran code, we overload the most often used and most numerically intense subroutines by linking them to BLAS and LAPACK. For example, we overload the matrix multiplications A*B by the BLAS matrix multiplication subroutines xGEMM, where x designates the data type of the matrices. This is done efficiently because the FortranMatrix object stores matrix element using Fortran convention. To make our code easy to use and understand, we name commonly used functions similar to Matlab commands. For example, to find the inverse of a matrix, we use B=inv(A), which is linked to the LAPACK subroutines xGETRF (LUfactorization) followed by xGETRI (inverse of an LU-factored matrix). Since our code spends most of the time doing matrix operations implemented with Fortran/Assembler native binaries, the overall performance of our code is close to a pure Fortran code. On the other hand, the codes are easy to write and understand. 3.4.3 BlockMatrix class We design a BlockMatrix class to implement mathematical algorithms as directly as possible. Without such an object class, block matrix algorithms are difficult to code and debug. In the following, we explain which features we want in a BlockMatrix class, and how these features enable us to implement complicated block matrix equations without much effort. We want to design a BlockMatrix class as a matrix FortranMatrix objects. Using the Fortran-Matrix object, the elements (blocks) of our BlockMatrix class will be intuitive and fast to operate. The storage should be efficient, which means that block matrices of different sparsity have different storage arrangements and different index functions. To do this, we define several subclasses such as the FullBlockMatrix, PeriodicBand-BlockMatrix, and the SparseBlockMatrix. The index functions are defined within the subclasses to overwrite the virtual index functions defined in the parent class. The object should be intuitive to use. For example, one may obtain the current operator from a tight-binding Hamiltonian with the following code: for (nt r=1; r<=N-1; r++) { J(r+1,r)=i*H(r+1,r); J(r,r+1)=dagger(J(r+1,r)); } Here, H(j,j+1) is a FortranMatrix itself, and also an element of the BlockMatrix H. In addition, we can access elements of the matrix H(j, j+1) by an expression like H(j,j+1) (ml,m2). These features allow us to implement block matrix algorithms efficiently and more importantly, debug in a shorter time. In the following, we show how the of Godfrin inverse method is implemented using the BlockMatrix class. We want to find G = M1, where M is a tridiagonal block matrix. These matrices are stored as PeriodicBandBlockMatrix, which is a subclass of BlockMatrix that stores a banded block matrix with periodic boundary conditions. For example, M is declared by 66 PeriodicBand-BlockMatrix< FortranMatrix< complex< double > > > M(N, KL, KD, KU, p, q); where N is the number of blocks. KL = KD = KU = 1 are the number of bands in the lower triangle, diagonal, and upper triangle or the matrix, and the block size is p x q. The Godfrin method shown in Equations A.21 and A.11 finds the inverse of a tridiagonal block matrix M in O(N2). To implement this algorithm, we first declare the subroutine, with input matrix M, and output matrix G and other variables used as temporary memory. Values such as the matrix size and block size are initialized, as shown below.3 template Godfrin(const PeriodicBandBlockMatrix< T > &M, PeriodicBandBlockMatrix< T > &G, FullBlockVector< T > &X, Full-BlockVector< T > &Y, Full-BlockVector< T > &C, FullBlockVector< T > &D) { int N=M.num-colsO; int p=M.blockho; int q=M.block-wo; T W(p,q); Here, class T is the class of an element of a block matrix, which is FortranMatrix< complex< double > > in this example. 3For clarity, the bound checking codes are not shown. Also, the actual code is slightly modified to output either the full matrix or selected bands of G depending on the input. 67- The next code segment is directly translated from Equation A.21. X(N)=0.0; for (int j=N; j>=2; j--) { W=inv(X(j)-M(j ,j)); C(j)=W*M(j ,j-1); X(j-1)=neg(M(j-1,j)*C(j)); } Y(1)=0.0; for (int j=1; j<=N-1; j++) W=inv(Y(j)-M(j,j)); D(j)=W*M(j,j+1); Y(j+1)=neg(M(j+1,j)*D(j)); Since in Equation A.21, Wj ans Zj are not used later, the memory of W is reused. Notice that variables like M(j-1, j) and W are themselves matrices. Although M is tridiagonal, there is no need to worry about where M(j-1, j) is stored when writing this code. The details of the storage are handled in the index function of the BlockMatrix object class according to the sparsity of the M. What is left is straight forward translation of equations. This code segment should be transparent even to people without C++ experience, and it is therefore easy to debug. 68 The last part of the subroutine computes the Green's function using Equation A.11. Again, this is a direct translation. for (int j=1; j<=N; j++) { G(j,j)=inv(M(j,j)-X(j)-Y(j)); for (int jl=j+1; jl<=N; jl++) G(j1,j)=C(j1)*G(j1-1,j); for (int jl=j-1; jl>=1; j1--) G(j1,j)=D(j1)*G(j1+1,j); } } The time spend on designing the BlockMatrix class should be less then the time we saved in coding and debugging other part of the code. Since the BlockMatrix object classes developed in this project are reusable, future implementations of mathematical algorithms involving block matrices should be as easy as the example shown above. CHAPTER 4 APPLICATION TO GIANT MAGNETORESISTANCE In the previous chapter, we explained how to compute the transport properties in magnetic nanostructures. In this chapter, we apply the method to compute the current-perpendicular-to-plane giant magnetoresistance (CPP-GMR) [4-14] in magnetic multilayers. As explained before, a good theory has to agree with both the conductivity and the GMR. To date, we have found two sets of experiments on the Fe/Cr CPP-GMR with both published resistivity and GMR data, one is by Gijs et al. [4] and the other is by Cyrille et al. [6-7]. In the following, we first review these experiments. The we try to see if we can use theoretical models to reproduce the experimental data for reasonable values of the parameters. In order to understand the trends shown in the experiments, the magnitude of the computed resistivity be similar to the values found in the experiments. The changes in resistivity and GMR as a function of different parameters are much smaller in scale, they should be studied after the magnitude of the resistivity is understood. Using the model in Chapter 2, we are able to calculate both the resistivity and GMR values which are close to Cyrille's data within reasonable parameters. However, the resistivity data of Gijs seem to be too small to fit in our theory unless very small disorder is used. We believe the required disorder is too small to be physical. 69 70 4.1 Experiment of Gijs and Collaborators The Fe/Cr superlattice in the experiment Gijs et al. [4-5] are micro-structured pillars, shown schematically in Figure 1.2. The pillar in the middle in Figure 1.2(d) is the sample to be measured. The samples sizes and measured conductivity and GMR are summarized in Table 1. Table 4.1: Summary of Gijs [Fe(30 A)/Cr(tcr)ioo samples Cr thickness tCr Area GMR psat(T = 0) A pm2 p2cm 10 90 1.08 11.0 28 20 0.08 24.3 40 6 0.04 25.2 The data are presented in terms of GMR = Ro/Rsat - 1, shown in Figure 4.1, and the quantity (Ro - Rsat)RoA2, shown in 4.1. Using these two quantities and the areas given in table 4.1, we can obtain the resistance Ro and Rsat., and hence Po and Psat. At low temperature, the resistivity of the samples is about 11 pA2cm for Cr layer thickness tcr = 1OA, and is increased to about 25 pf2cm for tcr = 28A and 40 A. The GMR of the sample is 108% for tcr = 10A and drops to 4-8% for tCr = 28A and 40 A. 71 As we shall see, the resistivity for both saturation field and zero field are very important to comparing with the theory. Iu I 100 80 60 40 20 i I 4 3 2 0 0 100 200 300 T(K) Figure 4.1: Temperature dependence of CPP and CIP GMR of pillar structures by Gijs et al. [4]. 4.2 Experiment of Cyrille and Collaborators Cyrille et al. [6-7] studied the effect of the CPP-GMR as the number of bilayers N of the Fe/Cr superlattice increases. They found that the GMR increases as N increases. 0 1 1 . 8 100 x (3 nm Fe + 2.8 nm Cr) SCPP (b) 4 -CIP 2 0 - 100 x (3 nm Fe + 4 nm Cr) (C) CPP -CIP 100 x (3 nm Fe + 1 nm CrO - (a) -I tc,(nm) -2.8 a P0.1 E 0.01 60 cc 40 20 0 0 1 00 200 30 T(K) Figure 4.2: Temperature dependence (Ro - Rsat)RoA2 for in the data of Gijs et al. [4], where A is the area of the sample. This quantity is used in the Valet-Fert theory [16]. Figures 4.2 show the resistivity and the GMR of a series of [Fe(30 A)/Cr(12 A)]N superlattices, where N ranges from 10 to 60. In this figure. the resistivity for both saturation field, psat, and zero field, po, is plotted as a function of the number of bilayers N. The magnitude of the CPP resistivity is about 35 PQcm for for saturation field, and is about 40-50 pQcm for the zero field. The experiments also show that Psat is roughly independent of N and po increases as N increases. This trend is not very clear on Figure 4.2(a), probably due to sample variations. However, since po and Psat are correlated, the trend is more 73 50- (a) 45( S40 35 30 0. 25 I I 40S 30- l 201010 20 30 40 50 60 N Figure 4.3: CPP Resistivities and GMR in FeCr superlattice as a function of number of bilayers measured by Cyrille et al. [6-7]. In (a), the zero field resistivity po increases as the number of bilayers N increases while the saturation field resistivity Psat remains roughly constant. In (b), the CPP-GMR increases from about 25% to 40% when N increases from about 17 to 40. This is a 60% change in the GMR. 74 clearly seen from the difference in resistivities which shows up in the GMR data in 4.2(b). If one concludes from the Figure 4.2(a) that psat is roughly constant, then po must be an increasing function of N. By detail analysis on the property of the samples, the authors show that the long length scale (about 10 nm) roughness of the interfaces increases as N increases, and suggest that the increase in po but not psat is due to the increase in spin dependent scattering at the interface. Since we only include atomic scale interface roughness in our calculation, we would not give the same trend as the experiments. 4.3 Comparing Theory with Experiments In Figure 4.4, we have plotted both the data from Gjis et al. experiment and those from the Cyrille et al. experiments on the same graph. Cyrille's data appears to be on the y-axis because they are measured at low temperatures, which is required for the leads to be superconducting. As we can see, the Gjis' resistivity is about a factor of two larger than the corresponding data of Cyrille. Our calculations are at zero temperature so we do not explicitly have temperature effects; however, within our model, we have several disorder parameters. By increasing the bulk disorder parameter, we can increase scattering, which is one of the effects of raising the temperature. In Figure 4.3, we have plotted the the maximum and saturation resistivity for our calculations as the spin-independent bulk disorder varies. A moderate residual bulk disorder is also added to the system. The residual bulk disorder is chosen so that resistivity calculations of bulk Fe and Cr films agree with experiments. The results in Figure 4.3 are close in magnitude to Cyrille's data and about twice larger than Gijs' data. 75 60 50 E40 0L0 Cn 20 . 10* 100 200 300 400 500 600 T(K) Figure 4.4: Resistivity in Cyrille's [Fe(30 A)/Cr(12 A)]N samples and Gijs' [Fe(30 A)/Cr(10 A)]ioo sample plotted on the same graph. The saturation field (triangles) and zero field (squares) resistivity in Cyrille's experiment are close the the y-axis because the experiments is done at the liquid helium temperature. The saturation field (pluses) and zero field (circles) resistivity in Gijs' sample are about 50% smaller than those in Cyrille's experiments. In order to fit Gijs' results, we have to calculate the resistivities using zero residual bulk disorder. The resistivities calculated this way are close to Gijs' data. However, the assumption that the bulk of the sample is perfectly grown is unreasonable. There are other disorder parameters related to the interface that we can change. In this calculation, we assumed that the interface thickness to be 4 monolayers, and the interface mixing percentage to be 20%. The resistivities can be reduced 7 6 60 50 e S U E40 Cn 30 201 0 100 200 300 400 500 600 spin-independent bulk disorder x 4 x 10 4 Figure 4.5: Calculated resistivity in [Fe(28.5 A)/Cr(11.5 A)],,, for moderate residual bulk disorder. Both the resistivity for parallel (squares) and anti-parallel (crosses) increase with the bulk disorder, however, the parallel increases faster, therefore the GMIR will decrease with the bulk disorder. The theory data is close in magnitude to Cyrille's data and about twice larger than Gijs' data. by reducing the interface disorder. However, we doubt that the interface can be much sharper than 4 monolayers and have less interface mixing. A better way to look at the differences among the data is to plot the change in resistivity against the saturation resistivity, as shown in Figure 4.3. For a fixed temperature and sample, there are two independent data, the resistivity at zero field and saturation field. By plotting functions of these two independent variables against each other, we can identify trends and remove the uncertainty in temperature dependence in the calculation. As seen in Figure 4.3, Gijs' data and 77 the theory have the same trend, due to a change in temperature or a change in bulk disorder, Cyrille's data represents the change in number of bilayers and has a different trend. 201- E ~15 =L x10 E 51- 01 0 Gijs a Theory Schuller a Xa a K x a X 0 Xa 3t 10 20 3 PsatA M 40 50 60 Figure 4.6: Change in resistivity against saturation resistivity. The data of Gijs (crosses) and the theory (squares) have the same trend; while the data of Cyrille (circles) has a different trend because it is done by varying the number of bilayers. Another advantage of plotting the data this way is when we compare two curves on the graph, we are comparing parameters other than the bulk disorders, which are the most difficult to model because they are usually not directly measured in the experiments. With the most uncertain parameter eliminated, a good theory should fit the experimental curve well. I As seen in figures 4.4 and 4.3, Gijs' data has a lower resistivity than both the theory and the experiment. The trends in Figure 4.3 do look similar, but there is a large offset between the curves. Fred Sharifi suggested one possibility is that there may be other conduction paths that allow electrons to by-pass the Fe/Cr superlattice. In Figure 4.3 we can obtain Gijs' curve by adding a shunt resistor to the calculated superlattice model. Superlattice shunt resistor Figure 4.7: Parallel resistor model. The effective resistivity is found by 1/Peff 1/piatt + 1/Pshunt, where Platt and Pshunt are the resistivity of the superlattice and the shunt resistor, respectively. In Figure 4.3, we plotted the the change in resistivity, po - psat, against the saturation resistivity for both Gijs' data and our calculation with and without the shunt resistor. We found that the calculation with a temperature-independent, 80 pQcm shunt resistor agrees well with Gijs' result. Therefore, Gijs data can be explained if in the experiment, some electrons are allowed to by-pass the superlattice and thereby reducing the measured resistivity. 79 18 1614. E 12 C.) _=10 CL8E 6. 4 2 0 I- 0 10 20 30 40 50 Psat (Q cm) Figure 4.8: Change in resistivity, po - Psat, against the saturation resistivity, Psat. The theory data without the shunt resistor (circles) are different from that of Gijs [Fe(30 A)/Cr(10 A)] (crosses), however, the theory data with the shunt resistor (squares) agrees with Gijs results very well. In order to fit the data, the resistivity of the additional path has to be almost independent of temperature and non-magnetic. It is not clear what caused the extra path. One possibility is that pin-holes are formed in the insulating polyimide, which allow conduction paths of Au to be deposited parallel to the sample. After the polyimide are washed away, the Au conduction paths are left on the sample. Another possibility is that the conduction paths may come from impurities at the edge of the superlattice. In Figure 4.9, we compare the theory with Gijs' [Fe(30A)/Cr(28A)] data. Again, with the shunt resistor, the theory agrees with the experiment very well. N x Gijs a Theory : shunt Theory: no shunt a 'S NU 1 N a N a a The shunt used in this case has a resistivity of 75 pAcm, this value is very close to the one obtained for the Fe(30 A)/Cr(1O A) sample. This suggests that the two paths are made of similar materials. By using our superlattice model with an extra conduction path way, we can explain both resistivity and GMR data for two different samples in Gijs' experiment. 1 1 1 E a- 8 6 4 2 0 S 10 20 30 Psat (gQ cm) 40 50 Figure 4.9: Change in resistivity, po - psat, against the saturation resistivity, psat. The theory data without the shunt resistor (circles) are different from that of Gijs [Fe(30 A)/Cr(28 A)] data. (crosses), however, the theory data with the shunt resistor (squares) agrees with Gijs results very well. Gijs - Theory: shunt Theory: no shunt 8 6 4 0 --'--------'-- 81 4.4 Other Effects on the GMR In this Chapter, we have assumed the adjacent magnetic layers to be completely anti-parallel at zero magnetic field. In experiments, the adjacent magnetic layers may not be completely anti-parallel due to the interlayer coupling [58-62]. The interlayer coupling controls the angle between the orientation of magnetization of adjacent layers at zero field. In Figure 4.10, the experiment of Parkin shows the saturation field oscillates as the thickness of the non-magnetic spacer layer varies. This effect also means that the magnetic coupling oscillates as a function of the layer thickness. The origin of the oscillation is the RKKY [58-59, 63-64] interaction. The GMR is also shown to oscillate with the thickness in the experiment of the Michigan State University group [9-10]. The effect is due to the GMR affected by how the magnetic layers align. 2.0 1.0 CO/V 0.5- 1.2 - (b) 0 4 .8 30. 10.4 10.0 - (c)7.5 - Co/Rh 5.0 _ -oR 2.5 -0 10 20 30 Spacer Layer Thickness (A) Figure 4.10: Dependence of saturation field on spacer-layer thickness for families of Co/V, Co/Mo, and Co/Rh multilayers [59]. 82 F200 e(5nm) (Co(1.5nm) /Cu(d) o CPP (HO) X CPP (HM) a CIP (HM) 150 100 50 CC 0x 0.0 2.5 5.0 7.5 10.0 dCu (nm) Figure 4.11: The magnetoresistance of Co/Cu layers as a function of the thickness of the space layer (Cu) for current perpendicular to the plane of layers (CPP) [9- 10]. Another complication is due to the fact that the magnetic layers are adjacent to non-magnetic layers or the vacuum. In such case, the magnetization in both the magnetic layer and the non-magnetic layer can be affected. More detailed information can be reviewed by studying the spin-density, such as those by Ohnishi et al. [42]. For example, Figure 4.12 shows that the spin density on the Fe layers extends to the vacuum, and there are negative spin density regions. However, since the computational cost in this effect is high, the GMR is usually computed 83 by assuming the magnetic layers are perfectly aligned in the saturation field and perfectly anti-aligned in the zero field. SPIN DENSITY Fe (001) S S- C Figure 4.12: Self-consistent spin-density map of seven-layer Fe(001) in unit of 0.0001 a.u. on the (110) plane. Each contour line differs by a factor of 2. Dashed lines indicate negative spin density [42]. CHAPTER 5 OPTIMIZATION OF GMR STRUCTURES In this chapter, we discuss the search for the optimal atomic configuration for the giant magnetoresistance. In previous chapters, we have discussed how to model the transport properties of magnetic nanostructures and compared the theory to experiments. Given a good model, we should be able to make predictions as well. In particular, what atomic configuration gives the highest GMR? This question is difficult because the configuration space is large. Take magnetic multi-layers made from 3d transition metals as an example, suppose we were to make a multilayer that consist 50 monolayers, and we can choose to use among Co, Ni, or Cu atoms in each monolayer, there are roughly 7 x 1023 possible configurations. It is hard to tell which configuration gives the highest GMR. Although theories are capable of computing the GMR for any of the configuration, it is impossible to exhaust every possibility. The conventional way of solving such a problem relies on physical insights, innovative ideas, and often trial-and-error combined with good luck. For a complicated function in the large configuration space, such as the GMR, explore within the subspace is limited to one that is simple enough to be understood with physical insights. In recent years, the combinatorics [65] method has been a useful tool in searching for new systems for interesting physical properties. For example, by systematically studying samples on a two-dimensional grid of gradual gradient of compositions, one can search though the configuration space at a much faster pace than was previously possible. However, trial-and-error 84 85 methods are unlikely to locate solutions in, for example, the 50-dimensional space mentioned. In order to search for the solution in a high dimensional space, the computer program must not search the configuration space randomly, but also make use of the local information. Optimization methods such as simulated annealing and genetic algorithms can be used. For example, using simulated annealing, Franceschetti and Zunger [37] are able to find the largest bandgap configurations in Al0.25Gao.75As alloys in a unit cell with 128 atoms by searching for about 104 steps, which is much smaller than the size of the configuration space. As discussed in Chapters 2 and 3, the GMR is complicated function of the atomic configuration and the disorder parameters. In the following chapter, we first review the simulated annealing method. Then we will apply the method to find out which configuration has the highest GMR, and discuss the finding and its implication. It is hoped that experimentally one can grow samples according to the optimal GMR configuration with non-equilibrium growth techniques such as the molecular beam epitaxy (MBE). 5.1 Inverse Problems Let f be a function of the atomic configuration C. We may want to find the minimum (or maximum) of the function in the configuration space fmi. = min f (C), C (5.1) 86 and the corresponding configuration Ci. In other situations, we may want to find the configuration that satisfies the target condition f(C) = ftrget, which can be found be minimizing of the absolute value of difference. More complicated requirements, such as prescribing multiple physical property, can also be formulated as the minimization of the function [37] 0(C) = Iw3 f (C) - flarget, (5.2) where wj is the weight assigned to the distance of property f, and the target target This problem belongs to a large class of "inverse problems". These problems are difficult because the functions in a high dimensional configuration space are complicated, and intuitive pictures are often not available. Finding the highest GMR configuration is an example of an inverse problems. In order to solve the inverse problem numerically, one must evaluate the function, or solve the "direct problem," many times. In practice, the direct problem itself is often difficult to solve. For example, it costs about 5 x 10" floating point operations to find the GMR of a superlattice with a unit cell of 50 atoms. To solve inverse problems, numerical optimization algorithms such as simulated annealing [32] and genetic algorithms [33-34] are used. These optimization techniques have been used to solve physics problems such as finding the equilibrium structure by minimizing the energy. However, the optimization describe in this chapter is rather different in nature. In solving the inverse problem, the solution obtained satisfies a set of conditions we prescribe, not the condition set by nature. This allow us to tailor-make a configuration with the desired properties rather 87 than settle for what is naturally available. The optimized configuration may be realized by growth techniques that can build structures monolayer by monolayer, or even atom by atom. In addition, the configuration space is very different from those used in the finding equilibrium structures. Here, we can change the chemical composition in the configuration, which is not a thermodynamical process. 5.2 Simulated Annealing In this section, we explain how to maximize a function by the simulated annealing algorithm. The flow chart of the algorithm is shown in Figure 5.1. To maximize a function f as a function of the configuration C, an initial configuration C is randomly generated and an initial "annealing temperature" T is set. The value of f(C) is then evaluated. A new configuration Ctrial is then generated by a Monte Carlo move from the previous configuration, and the value of f(Ctriai) is also evaluated. If the functional value of f for the new configuration is higher than that for the previous configuration, it is accepted. On the other hand, if the new GMR is lower, it is only accepted with a probability of exp((f(Ctria) - f(C))/T). The process is continued as the annealing temperature is gradually lowered until there is no further change in the configuration and the final configuration is a near-global maximum. This process resembles the annealing process in statistical mechanics, where the total energy is minimized when a system is cooled slowly. In simulated annealing, the value of f is used to decide whether to keep or discard the trial configuration. Compared with random walks routines, simulated annealing causes the system to spend more steps near the configurations possessing large values of f. Simulated annealing is easy to implement, and it is typical that 88 the number of configuration tested in simulated annealing be much smaller than the size of the configuration space. Set Initial configuration C and pseudo temperature T Monte Carlo move Cto Cena FInd Af =RuCtmiaet -f (C) No Af > 0? Yes Accept Cimia with a probability Accept CGnal Pe*A Slowly decrease T to 0 e.g. T =Toexp(-'nrj) No -O Ye s T ~N o? E n d Figure 5.1: Flow Chart for Simulated Annealing. The flow chart for using simulated annealing to search for the optimal configuration C, that maximize a function fis shown. 5.3 Accuracy and Speed Requirements on Model Searching for the optimal structure using optimization techniques requires a delicate balance between accuracy and speed. On one hand, the model has to be accurate enough so that the optimal solution is reliable. The model has to capture the physics and to accommodate a wide range of configurations. The level of detail in the model varies from problem to problem. For optimization of the transport properties for nanostructures, the model has to include important effects of the band structure of the materials, the geometry of the nanostructures, and the scattering effects caused by impurity atoms and structural disorders. All of these factors significantly affect electron transport. This requirement eliminates methods containing less microscopic detail than the multi-band tight-binding model that we discussed in previous chapters. On the other hand, calculations of the transport properties in this model have to be fast enough so that it can finish in reasonable time. The number of annealing steps is expected to be of the order of 103 to 104. To finish the optimization within a month, one has to be able to compute the function for a configuration in about 4 to 40 minutes. This is very demanding even for running tight-binding codes on a computer cluster. At first glance, the matrix operations of finding the Green's function, vertex corrections, and conductivity in Chapter 2 are are all O(N3) operations, where N is the number of atoms in the unit cell. In addition, the self-energy and the vertex corrections are solved by self-consistent iterations. Performing optimization with a self-consistent O(N3) calculation would take about 50 CPU-years, which is impractical. This is why we have developed the fast algorithm in Chapter 3. For N = 50, this algorithm is more than 100 times faster than the O(N3) algorithm. Running on a Linux cluster with ten 866MHz Pentium III CPUs, it takes about 20 minutes for each GMR calculation and about two weeks for an optimization with 1000 annealing steps. 90 The computational cost counting in Chapter 3 are useful in estimating the computer resources required. A fair estimation of the time for an optimization is required before coding could be started. 5.4 Constrained Optimization of the GMR In this section, we explain how to apply simulated annealing optimizing the giant magnetoresistance. As discussed in Chapter 3, there are two definition of the giant magnetoresistance. Here we adopt the "inflated" definition of the GMR: GMR - 1, (5.3) O'AP because it is "optimization friendly". If we were to use the alternate definition, GMR = 1 - U'AP/OP, the GMR would be bounded by 100% and a significant improvement from, say, 80% and 85% would look small compared to the annealing temperature, and could be washed away until the temperature is very low, hence extending the annealing time significantly. When the physical property to be optimized is "optimization unfriendly," one needs to find a monotonic function of the property that magnifies the regime of interest. Not every configuration qualifies as a GMR structure. Unqualified configurations, such as those with no magnetic layers, will have zero or very small GMR. When the computational cost for each configuration is small, one can rely on the exponentially small probability to make sure the unqualified configurations are not chosen. However, when the computational cost for each configuration is large, the time wasted in computing the unqualified structures become very costly. In com- 91 puting the transport properties, self-consistencies in the self-energy and the vertex correction are required, therefore the result for the last configuration is not useful in finding the result for the present configuration. Therefore each GMR calculation is expensive. One can also use a test algorithm to reject unqualified configurations, but it would still waste significant amount of time when the qualifying configuration subspace is much smaller than the full configuration space. Therefore, it is better to perform the optimization in the constrained subspace of GMR structures. The algorithm of generating the configurations in the constrained subspace will be more complicated, but it will produce results in a much shorter time. To optimize the GMR using simulated annealing, we consider the constrained subspace of configurations consisting of a unit cell with two magnetic layers separated by two non-magnetic spacer layers. After an arbitrarily chosen initial unit cell, each subsequent configuration is generated by one of the following Monte Carlo moves: (i) inserting a monolayer, (ii) removing a monolayer, or (iii) changing the composition of a monolayer. To make sure we do not leave the constrained subspace of GMR, a move is rejected unless it will not eliminate a layer, or turn a magnetic layer into a non-magnetic layer by, for example, changing the atoms from Co to Cu. To prevent complications such as pin holes and magnetic dead layers, we further limit the minimum thickness of a layer to two mono-atomic layers. If the GMR of the new configuration is higher than the previous configuration., then it is accepted. On the other hand, if the new GMR is lower, it is only accepted with a probability of exp(AGMR/T), where AGMR is the change in the GMR and T is the simulated annealing temperature. The process is continued as 92 the annealing temperature is gradually lowered until there is no further change in the configuration and a final near-global maximum in the GMR is attained. 5.5 Optimal GMR Configuration in CoNiCu superlattices In the following, we discuss the parameter and results of a study looking for the optimal configuration for superlattices made with Co, Ni, and Cu, in conditions similar to GMR experiments at the room temperature. To make the study more manageable, we only consider the fcc(111) lattices. As mentioned in Chapter 2, the tight-binding parameters for Equation 2.1 are obtained from fits to density-function calculations where up to second-nearestneighbor hoping energies are included. In order to simulate experimental conditions, we fixed the bulk, spin-independent scattering parameter Vbulk at 0.2 eV2. Using this parameter, the resistivity for bulk Cu is calculated as 3.1 [Mcm, about twice that of the room temperature resistivity of Cu. Therefore, we are modeling a disordered metal system. The spin-dependent impurity potential caused by interface diffusion is found by fitting the tight-binding model to a supercell density-functional calculations of an impurity atom in the host. The values for the impurity potentials obtained are as follows: for the majority spin, Auco/cu =1.7 eV, AUNi/Cu =0.56 eV, AUC0/Ni =0.56 eV; and for the minority spin, Auco/Cu =2.9 eV, AUNi/Cu =1.05 eV, AUCo/Ni =0.64 eV. The interface mixing percentage is fixed in the follow way: for each monolayer, 10% of the atoms diffuse into each of the two monolayers above and below. Using the above parameters, the GMR for [Co3/Cu3/Co3/Cu3] superlattice is 36%. This value is similar to experiments. The exact value for the mixing percentage in the experiment is unknown because it is |

Full Text |

PAGE 1 ELECTRON TRANSPORT THEORY IN MAGNETIC NANOSTRUCTURES BY TAT-SANG CHOY A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2001 PAGE 2 ACKNOWLEDGMENTS I would like to thank Selman Hershfield for being my thesis advisor. He introduced me to the field of magnetic nanostructures and taught me Green's function techniques. I also learned from Selman how to work with experimentahsts. He has given me guidance in research and also encouraged me to find research topics myself. When we started working on optimizing the GMR structure, I was not optimistic that the difficulties could be overcome. Without Selman's encouragement and support, I would have waited a few years before taking this challenge. I would like to thank Professors Arthur Hebard, Peter Herschfeld, Mark Meisel, Stephen Pearton, Fred Sharifi, and Christopher Stanton for serving as as members of my supervisory committee. Fred Sharifi taught me a lot about sample growth, information I used in my research. I also enjoyed a lot the ad hoc discussion group organized by Professors Arthur Hebard, Selman Hershfield, Mark Meisel, Andrew Rinzler, and Fred Sharifi. Two good friends influenced me greatly in my physics career. Hoi-Fung Chau was the first research physicist I have known. Ever since high school, I was lucky to have learned from him many new ideas in physics, mathematics, and computer sciences. Those ideas helped me realize my true passion is in physics rather than computer science. Jian Chen taught me a lot during and after he worked in Selman's group as a postdoc. I enjoyed dropping by his office chatting with him about physics and other things. I regard Jian Chen as my first close collaborator in science. I would like to thank fellow graduate students Kingshuk Majimidar, Steven Patamia, Vladimir Boychev, Zhihong Chen, Stephanie Getty, Heather Hudspeth, ii PAGE 3 Quentin Hudspeth, and Kelvin McCarthy for their friendship. Steven Patamia gave me useful advice in C++ programming. Many thank to our former system administrator Chandra Chegireddy, who helped tremendoiisly in setting up the special hardware and software for my research. I would like to thank Dr. Michael Jones for making the LaTex dissertation template available. The work in this thesis was supported by grants from the Air Force Office of Scientific Research, the National Science Foundation, and the Research Corporation Company. Finally, I would like to thank my parents and my sister for their support. I am grateful for my parents' understanding of their son spending the past few years eleven time zones away. Also, my sister's gifts from time to time are much appreciated. iii PAGE 4 TABLE OF CONTENTS ACKNOWLEDGMENTS ii ABSTRACT . vii CHAPTERS 1 INTRODUCTION 1 1.1 OMR Experiments 2 1.2 Models of the CPP-GMR 4 1.2.1 Resistors-in-Series Model 5 1.2.2 Fermi Surface Mismatch Picture 7 1.2.3 Potential Step Picture 10 1.2.4 Effects of the Disorder 12 1.3 Optimization at the Atomic Level 14 1.4 Overview of the Rest of the Dissertation 17 2 MODELING TRANSPORT IN NANOSTRUCTURES 18 2.1 Modeling of the Band Structure 20 2.1.1 Multi-band Tight-binding Hamiltonian 20 2.1.2 Modeling Impurities 21 2.1.3 Impurity Averaged Green's Functions 23 2.2 Conductivity Calculation 26 2.2.1 Formal Aspects of Cmrrent Operators 27 2.2.2 Cmrrent Operators in One Dimension 31 2.2.3 Kubo Formula 36 2.2.4 Cm-rent Conservation 39 2.3 Realistic Parameters 43 2.3.1 Tight-Binding Parameters 43 2.3.2 Imptuity Potential for Substitutional Disorder 44 2.3.3 Impurity Potential for Structural Disorder 46 2.3.4 Extension to Magnetic Alloys 46 iv PAGE 5 3 COMPUTATION PROCEDURE 49 3.1 Current-Perpendicular-toPlane Conductivity 50 3.1.1 Impurity Averaged Green's Function 51 3.1.2 Vertex Correction 53 3.1.3 Conductivity 54 3.1.4 Diagonal Impurity Potential 56 3.1.5 Total Computational Cost 58 3.2 Current-In-Plane Conductivity 59 3.3 Conductivity Tensor 60 3.4 Coding ,^ 61 3.4.1 Template Numerical Toolkit (TNT) ... 62 3.4.2 Fortran performance with BLAS and LAPACK 64 3.4.3 BlockMatrix class '.' 64 4 APPLICATION TO GIANT MAGNETORESISTANCE 69 4.1 Experiment of Gijs and Collaborators 70 4.2 Experiment of Cyrille and Collaborators 71 4.3 Comparing Theory with Experiments 74 4.4 Other Effects on the GMR 81 5 OPTIMIZATION OF GMR STRUCTURES 84 5.1 Inverse Problems 85 5.2 Simulated Annealing 87 5.3 Accuracy and Speed Requirements on Model 88 5.4 Constrained Optimization of the GMR 90 5.5 Optimal GMR Configuration in CoNiCu superlattices 92 5.6 Understanding the Large GMR 94 5.7 Conclusion 97 6 CONCLUSION 100 APPENDICES A FAST INVERSE ALGORITHMS 103 A.l Godfrin Inverse Algorithm for Tridiagonal Block Matrix 104 A.2 Inverse of Tridiagonal Block Matrix with Periodic Boundary Condition 108 V PAGE 6 A.2.1 Rank-one changes: Woodbury formula 108 A.2.2 Applying the Woodbury Formula 109 A. 2. 3 Cost Considerations 110 A.3 Efficient Coding for the Godfrin Algorithm 112 B COHERENT POTENTIAL APPROXIMATION 115 REFERENCES 118 BIOGRAPHICAL SKETCH 122 vi PAGE 7 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy ELECTRON TRANSPORT THEORY IN MAGNETIC NANOSTRUCTURES Tat-Sang Choy , . ' August 2001 Chairman: Selman P. Hershfield Major Department: Physics Magnetic nanostructure has been a new trend because of its application in making magnetic sensors, magnetic memories, and magnetic reading heads in hard disks drives. Although a variety of nanostructures have been realized in experiments in recent years by innovative sample growth techniques, the theoretical study of these devices remain a challenge. On one hand, atomic scale modeling is often required for studying the magnetic nanostructures; on the other, these structiures often have a dimension on the order of one micrometer, which makes the calculation numerically intensive. In this work, we have studied the electron transport theory in magnetic nanostructiures, with special attention to the giant magnetoresistance (GMR) structmre. We have developed a model that includes the details of the the band structmre and disorder, both of which are both important in obtaining the conductivity. We have also developed an efficient algorithm to compute the conductivity in magnetic nanostructures. The model and the algorithm are general and can be applied to complicated structinres. We have applied the theory to current-perpendicular-toplane GMR structures and the results agree with experiments. Finally, we have vii PAGE 8 searched for the atomic configuration with the highest GMR using the simulated annealing algorithm. This method is computationally intensive becaiise we have to compute the GMR for 10^ to 10^ configurations. However it is still very efficient because the number of steps it takes to find the maximum is much smaller than the number of all possible GMR structures. We found that ultra-thin NiCu superlattices have siurprisingly large GMR even at the moderate disorder in experiments. This finding may be useful in improving the GMR technology. viii PAGE 9 CHAPTER 1 INTRODUCTION Magnetic nanostructures [1-2] are magnetic devices which have at least one dimension of the order a nanometer, 10~^ m. These devices have apphcations in making magnetic sensors, magnetic memories, and magnetic reading heads in hard disks drives. Predicting the electrical transport properties theoretically is a challenge because the detailed arrangement of atoms affects the observable properties significantly. In this dissertation, we will develop methods to predict the transport properties of magnetic nanostructm-es by taking into account band structiue and disorder. We will first develop models to compute the giant magnetoresistance ratio (OMR), which is the fractional change of the resistance due to an applied magnetic field, of multilayer structures and compare our calculations with experiments. We will then solve the "inverse problem," which is to find the atomic configviration that has the highest GMR. Although the method we develop can be applied to various geometries, in this dissertation we concentrate on magnetic superlattices, which are multilayers formed by repeatedly stacking magnetic layers followed by non-magnetic layers. The thickness of each layer is of the order of one nanometer; however, the other two dimensions can be as large as one micrometer, 10"^ m. This geometry can be grown by techniques used in depositing thin films. In the following, we will first discuss the GMR experiments, explain several pictures for understanding the gi1 PAGE 10 2 ant magnetoresistance, and discuss the recent trend of optimization on the atomic level. 1.1 GMR Experiments ' The first "giant" magnetoresistance measurement was that of Baibich et al. [3], where they observe an almost 50% reduction in the resistance of Fe/Cr superlattices in an applied magnetic field of 20 kG. As shown in Figure 1.1, the resistance at zero field, Rq, is the highest. As the magnetic field is applied, the magnetization of the adjacent Fe layers align, which causes the resistance to decrease. At a field larger than the saturation field Hs, the magnetization of the adjacent Fe layers is parallel to each other. Therefore, the resistance stays at the value of the saturation field, Rsax. The giant magnetoresistance ratio is defined as GMR= (1.1) In this experiment, the current travels in the plane of the multilayers, which is referred to as the ciurrent-in-plane or CIP giant magnetoresistance. In this dissertation, we focus on another geometry, the current-perpendicularto-plane (CPP) geometry, in which the current fiows perpendicular to the plane of the multilayer. The CPP-GMR is usually larger than the CIP-GMR for the same lattice. Most experiments to date are for the CIP geometry, while a smaller portion of the experiments are for the CPP geometry [4-15]. Current-perpendicular-toplane experiments are difficult to perform because the thickness of a superlattice is of the order of a hundred nanometers at most, while the area of the layer is PAGE 11 3 Magnetic field (kG) Figure 1.1: Magnetoresistance of three Fe/Cr superlattice at 4.2 K [3]. The current and the apphed field are along the same [110] axis in the plane of the layers. Hg is the saturation field required to align the magnetizations of Fe layers. As the magnetizations are aligned, the resistance is reduced from the peak at H = 0. usually much larger than 1 //m x 1 //m. When a current passes though the system perpendicular to the plane, the resistance is too small compared with the resistance of the lead and the contact, and is thus hard to measure. To measiue the CPP resistivity, one must improve the aspect ratio of the sample. Small cross-sectional area to length ratio has been produced by using micro-structiured pillars (Figiire 1.2) [4], by growing superlattices inside the pores in polycarbonate (Figure 1.3) [15], and by connecting several superlattices with superconducting leads (Figure 1.4) [6-7]. PAGE 12 4 substrate ib) ' Mo 100 nm AI2O3 200 nm Mo 100 nm Au 300 nm Fe/Cr 500 nm Cr *00 nm Si02 (c) Polyimide Figure 1.2: Schematic diagram the samples growth process used by Gijs et al. [4]. The cross sectional area S ranges from 6 to 130 /im^. 1.2 Models of the CPP-GMR Â• There has been important progress in understanding the CPP-GMR effect. Classical resistors-in-series models can be applied to study the CPP-GMR when the layer thickness is much larger than the Fermi wavelength [16]. The origin of the spin anisotropy, which measures the difference between the resistivity in the two spins, is more difficult to understand. Various free-electron models [17-21], tightbinding models [22-26], and densityfunction methods [27] have been performed. In the following, we discuss several models for the CPP-GMR. In the CPPGMR system, each electron has to go through every layer and interface. It is thus easier to understand intuitively than CIP-GMR system, where electrons travels along the layers, while at the same time going across interfaces of layers and/or scattering from the interfaces. ^ PAGE 13 5 Figure 1.3: Schematic representation of the array of nanowires in the insulating polymer matrix [15]. 1.2.1 Resistors-in-Series Model When the layer thickness is smaller than the spin-diffusion length [28-29], one can treat the two spins independently, and this approach is called the two-fluid model [16]. When the layer thickness is much larger than the Fermi wavelength, the GMR effect can be explained with the classical resistors-in-series model, which assumes electrons pass through spin-dependent resistors, as shown in Figure 1.5. Let the resistance for the majority and minority spin in the bilayer formed by a magnetic layer and a non-magnetic layer be R and r, respectively. At saturation field, the magnetic moment of adjacent magnetic layers are parallel to each other. Assuming the spin diffusion length is much larger than the layer thickness, the resistance for the majority and minority spin is R + R and r + r, as shown in PAGE 14 Nb I I SiQ, (a) (b) Figure 1.4: CPP sample of Cyrille et al. [6]. (a) Schematic cross section of a CPP sample; (b) optical micrograph of a typical sample. Figure 1.5(a). Therefore, the net resistance for the saturation field is given by _ 1 Rr 2R + r (1.2) At zero field, the magnetic moment of adjacent magnetic layers are assumed to be anti-parallel to each other, as shown in Figm-e 1.5(b). A majority spin electron of the first magnetic layer will be in the minority spin channel of the second magnetic layer. Therefore, the resistance for each spin channel is r + i?. As a result, the net PAGE 15 7 resistance at zero field is given by i?o Â— R + r (1.3) Since (i? + r)^ > 4Rr, one can see that Rq is larger than i2sat, which is the explains the GMR. (a) (b) R R R R Figure 1.5: Simple resistor model for CPP-GMR. The circuits corresponding to (a) parallel, and (b) anti-parallel, magnetic configuration. In both cases, the majority and minority spin electrons pass through the left and right series of resistors, respectively. In (a), the resistance for the majority and minority spin channels are 2R and 2r, respectively; the net conductance is l/2R+l/2r. In (b), the resistance for either channels are r + R; the net conductance is 2/(r + 1.2.2 Fermi Surface Mismatch Picture Although the basic idea of the GMR can be imderstood by the simple resistor model above, this model does not explain the origin of the spin anisotropy, which is the difference divided by the sum of the resistance in the two spin channels. There are two spin anisotropics, the bulk spin anisotropy and the interface spin anisotropy. The interface spin anisotropy can be understood with simple pictures PAGE 16 8 and is often the dominant effect of the GMR unless the layers are thick. Therefore, in the following we will explain the GMR using only the interface spin anisotropy. The most intuitive way of explaining the spin anisotropies and the CPP-GMR is probably the Fermi surface picture. We use the Fermi surfaces plots from our "Fermi Surface Database" website at http://www.ufl.edu/fermisurface [SCSI] in the following discussion. It provides both pictures and interactive Virtual Reality Modeling Language (VRML) models of Fermi surfaces for most elemental metals. We assume the layers to be large compared with the Fermi wavelength such that the bulk Fermi surfaces is a reasonable approximation. Again the two fluid model is used, assuming the spin diffusion length is large compared to the thickness of the layers [16, 28-29]. As an example, we consider the Fe/Cr/Fe system, a three-layered structure formed by sandwiching a Cr layer in between two Fe layers. The Fermi surfaces of bulk Fe spin-up and spin-down channels, and non-magnetic Cr, are arranged in Figure 1.6 for both (a) the parallel and (b) the anti-parallel magnetic configurations. As shown in Figm-e 1.6(a), the Cr Fermi surface is very different from the Fe spin-up Fermi siuface, while it is much closer to the Fe spin-down Fermi surface. If an electron in the spin-up channel of the Fe layer on the left needs to go into the Cr layer, there will be larger scattering at the interface because a wavevector on the spin-up Fermi surface of bulk Fe is not on the bulk Cr Fermi surface. When this electron goes from the Cr layer to the Fe layer on the right, it is scattered again. Therefore, an electron goes though two large scattering events, and as a result, the spin-up channel has a low conductivity. On the other hand, a spin-down electron from the Fe layer on the left is not scattered as much in going PAGE 17 9 through the two interfaces. The scattering is much smaller because the Cr and Fe spin-down Fermi surfaces are similar. Thus, the spin-down channel has a higher conductivity. The total conductivity for the parallel configuration is the sum of the two channels, which is dominated by the spin-down channel. (a) (b) Spin Up Channel Cr Feup'T Spin Down Channel Fedn Cr Fedn Spin Up Channel -eup Cr Fe< t in Spin Down Channel Fedn Cr Figure 1.6: Fermi Siuface Mismatch in the Fe/Cr/Fe multilayer. The bulk Fermi for each spin in the layers in the (a) parallel, and (b) anti-parallel, magnetic configuration are shown. There are only three different Fermi surfaces, Fe majority spin, Fe minority spin, and Cr. Since the spin of an electron crossing the layers is assumed to fixed, the electron behave as a majority or minority electron depending on whether its spin is aligned or anti-aligned with the local magnetization of the layer. The anti-parallel case is shown in Figure 1.6(b). The spin-up and spin-down Fermi smfaces for the Fe layer on the right are exchanged. When the local magnetization on the left and the right are anti-parallel to each other, an electron which PAGE 18 10 starts as a spin-up (majority) electron on the left will be in the minority population on the right. Hence it behave as a spin-down electron in the bulk of the Fe layer on the right. Here, electrons in the both spin channel go though the same amoimt of scattering. Since the Cr/Fe spinup Fermi surface mismatch is much larger than the Cr/Fe spin-down Fermi siurface mismatch, the conductivity for the spin-down channel of the parallel configuration is much larger than the other channels. Therefore, the conductivity for the parallel magnetic configuration is much larger than the conductivity for the anti-parallel magnetic configuration. In this case, the giant magnetoresistance is positive. The CPP-GMR in Co/Cu/Co can be explained similarly using the Fermi surfaces in Figure 1.7. Here, the Cu Fermi surface is similar to the Co spinup Fermi surface, but rather different from the Co spin-down Fermi surface. Therefore, the conductivity is dominated by the spin-up channel of the parallel configuration. Again, the conductivity for the parallel magnetic configuration is significantly larger than the conductivity for the anti-parallel magnetic configuration. Hence, the OMR is also positive. We note that the spin anisotropy for the Fe/Cr interface has a different sign than the spin anisotropy for the Co/Cu interface. 1.2.3 Potential Step Picture A related way of understanding the spin anisotropy, or the difference in resistance for difference spins, at interfaces is to look at the potential step across the interfaces. As will be explained in Chapter 2, the major difference in the PAGE 19 11 (a) Co fccup J i Cu Cojccup Spin Up Channel [spin Down Channel Cu (b) Spin Up Channel Cojccup f Cu 1 jSpin Down Channel Cu \4f^/ Cojccup f { Figure 1.7: Fermi Surface Mismatching the Co/Cu/Co multilayer. The bulk Fermi for each spin in the layers in the (a) parallel, and (b) anti-parallel, magnetic configinration are shown. ' tight-binding parameters in the 3d transition metals is the on-site energy of the d-orbitals. The system looks like a simple transmission-though-potential-barrier problem. We must be careful, however, that the potential step affects only the d-electrons, while the sand p-electrons carry a significant fraction of the current. The effect is comphcated by the spd-hybridization. From the tight-binding parameters, the Cu potential is 0.82 eV lower than the spin-up Co potential and is 2.3 eV higher than spin-down Co potential. Since the potential step in Cu/Co spin-up interface is higher, it has higher resistance. PAGE 20 12 Similarly, the Cr potential is 1.53 eV higher than the spin-up Fe potential and is 0.23 eV lower than the spin-down Fe potential. The potential step may also be used to estimate the impurity potential due to interface diffusion of, say Cu into Co. ^ 1.2.4 Effects of the Disorder Although the simple pictures in the last section can help imderstand the origin of the giant magnetoresistance effect, these pictures do not include both the band structure effect and the disorder effect. We have explained that the band structure is often the source of the GMR. On the other hand, both the conductivity and the size of the GMR are determined by the disorder. If we compute the conductivity with small disorder, we would get a conductivity orders of magnitudes larger than what is usually obtained in experiments. Also, if the disorder is too small, the GMR would be too large compared to experiments. Therefore, disorder must be included in the calculation in order to obtain the conductivity and the GMR. As explained by J. Chen et al. [25], the disorder can affect the the density of states (DOS) significantly. The density of states of Cu and Co with different disorder parameters is shown in Figure 1.8. As we can see, the disorder significantly changes the density of states. The effect of the disorder is to spread out the DOS. The net effect at the Fermi surface is hard to predict without the calculations done in Chapter 2. In Cu and the spin-down channel of Co, the DOS at the Fermi level is not significantly affected. However, for the spin-up channel of Co, which has a similar band structure to Cu, the DOS is doubled due to the DOS spread from PAGE 21 13 the d-peak just below the Fermi level. We will see in the next chapter that this significantly reduces the conductivity, even though the DOS is doubled. -4.0 -2.0 0.0 E-Ep (eV) Figure 1.8: Density of states (DOS) of Cu, and Co spinup and spin-down channels for different impmity density. When the impurity density is small, the DOS remain close to the unperturbed DOS. However, when the impurity density is moderate, the DOS can change a lot. The conductivity is affected by the DOS near the Fermi level. For Cu and Co spin-down, the change in DOS at the Fermi level are small. However, for Co spin-up channels, which has a similar band structure to Cu, the DOS at the Fermi level is increased by a factor of two [25]. To calculate the GMR, we have to find the appropriate disorder parameters. We set the bulk impiurity parameter so that the bulk resistivity is similar to experiments. A single bulk impurity parameter should be sufficient to give the right bulk resistivity for all elements in a sample. In Chapter 1, we explained how to find PAGE 22 14 the impurity potential due to diffusion. Tlie density of the impurity is harder to set because there is no experiment the clearly measure this parameter. However, since only the total scattering rate matters, one can set the parameters to match the experiments. Since the same impurity density has to be used for both parallel and anti-parallel magnetic configurations, the single parameter should be able to give the correct resistivity for both configmration. In other words, there should less parameter than the niraiber of data points. In Figure 1.9, the current-perpendicular-to-plane resistivity times the layer thickness t for the Co(/Cut superlattices is plotted against t. The results for the spin-up channel, labeled as Cott/Cut, and the spin-down channel, labeled as Coit/Cuf, are shown for calculations with bulk disorder only, bulk plus spinindependent interface disorder, and bulk plus spin-dependent interface disorder. The lines are Uneax fits to the computed results. The interface resistances can be read from the intercepts. The bulk resistance is related to the slope. Both the interface resistance and the bulk resistance are affected by the types of disorder. It is found that both bulk disorder and spin-dependent interface disorder are necessary to explain experiment results [25] . 1.3 Optimization at the Atomic Level Optimizations in physical systems axe often performed by the two popular methods: simulated annealing [32] which is inspired by statistical mechanics, and genetic algorithms [33-34], which are inspired by evolution. Simulated annealing has been used in optimization problems such those in circuit design [32]. It is also used in problems such as structviral optimization. For example, it has been used PAGE 23 15 t(A) Figure 1.9: Current-perpendicular-to-plane resistivity of superlattices Co/Cu. The spinup and spin-down channels are labeled Cott/Cut and CoJ.t/Cut, respectively. The resistivity are calculated for three cases: with only bulk disorder, bulk plus spin-independent interface disorder, and bulk plus spin-dependent interface disorder. The lines are linear fits to the computed results. The interface resistances can be read from the intercepts [25]. to find the most stable structure of a 7 Si atom cluster using the pseudo-potential density-functional theory [35]. The genetic algorithm has also been used in optimization of the molecular structure. For example, Deaven and Ho [36] performed a structural optimization by starting with a population of random molecules with 60 caxbon atoms. Then the more stable, or "fit" , molecules are cut into half and glue, or "mate", with each other to produce a new population of molecules. Genetic manipulations such as crossing over are performed in the process. As the simulation goes on, the most stable structure, the ideal icosahedral buckyball structure, is achieved. In Chapter 5, we turn to a new direction of research in which we optimize the atomic configuration of a nanostructure to obtain the maximiun GMR. Performing optimization of physical properties in the atomic scale is an exciting new trend. PAGE 24 16 Pranceschetti and Zunger performed the simulated annealing optimization to find the atomic configurations with the maximum bandgap in semiconducting alloys and superlattices [37] . To maximize the bandgap energy with simulated annealing (Figure 1.10), they started with a random atomic configuration and performed Monte Carlo moves in the configuration space, which is the space of all possible configurations. A move can be accepted with a probability even when the new configuration has a smaller bandgap. However, as the annealing goes on, the probability is reduced until the system finally rests in the maximum bandgap configuration. The maximum bandgap configuration for a Alo.25Gao.75As alloy is shown in Figiu-e 1.11. It would be difficult to find this configuration using trialand-error. 1-90, , , , , , 1 1.74 1 1 1 ' 1 1 1 0 2.000 4,000 6,000 8.000 10,000 12,000 Monte Carlo move Figure 1.10: Optimization of the bandgap energy of Alo.25Gao.75As alloy by simulated annealing [37]. The bandgap energy fluctuates as the Monte Carlo routine progress. However, since there is a bias towards increasing bandgap and the amount of fluctuation is decreasing, the system finally arrive a maximum bandgap configuration. Although atomic level optimization has been performed, the optimization of physical properties other than the energy is new. Here one searches for an atomic PAGE 25 17 Figure 1.11: Maximum bandgap configuration of Alo.25Gao.75As [37]. This complicated configuration is unlikely to be found by trial-and-error. configuration that is, by definition, not the most stable configuration. However, with the non-equilibrium growth techniques, structures optimized for the desire physical property may become a reality. 1.4 Overview of the Rest of the Dissertation The rest of the dissertation goes as follows. In Chapter 2, we will develop a method to compute the conductivity in magnetic nanostructures. Both the band structure and the disorder will be included. In Chapter 3, we will discuss how to apply the model efficiently by rewriting the equations. In Chapter 4, we apply the method to calculate the CPP-GMR and compare the result with experiments. In Chapter 5, we will optimize the atomic configuration of a nanostructure to obtain the maximum OMR. In the final chapter, we will discuss the imphcations of this work and future directions. PAGE 26 CHAPTER 2 MODELING TRANSPORT IN NANOSTRUCTURES In this chapter, we discuss how to model electron transport in magnetic nanostructures. The common approaches include semi-classical methods, simplified one or two band tight-binding model, multi-band tight-binding model, and densityfunctional calculations. Each model contains an increasing level of details and has its advantages and disadvantages. Semi-classical methods [16] can be used for studying large systems with complex geometry that other methods can not solve efficiently. The conductivity is often calculated with the Boltzmann equation. Bulk and interface scattering effects can be included as parameters in the Boltzmann equation, without the need for microscopic details. This approach valid when the length of smallest feature is much larger than the Fermi wavelengths, but it is often possible to extend it to smaller systems for qualitative study. However, when the feature length scale is only a few nanometers, the semi-classical method may become problematic. The simplified one or two band tight-binding model is often used to give qualitative understanding of the finite size effects and some of the band structiure effects. Usually, the parameters in this model axe obtained by reasonable estimations. The conductivity is usually obtained by the Kubo formula [38] or the Landauer formula [39-40]. This model is easy to solve, and sometimes it is even possible to obtain analytical solutions. In addition, the simplified band structure is helpful in understanding some effects that may be less obvious in more complicated 18 PAGE 27 : .. 19 models. However, in real materials, such as iron or nickel, the band structxure is too complicated for the results of this simple model to compare with experiments quantitatively. The multi-band tight-binding model [22-26] is used for quantitative comparison with experiments. Parameters are often fit from density-functional theory calculations and hence it captures the band structure. Details of atomic variation, impurity scattering effects [25-26] can be modeled accurately. Although this method is not as accurate as a full density-function theory calculation, it has the advantage of being much faster. It is also slower than the simple tight-binding model; however, as we will show in detail in Chapter 3, it is fast enough to perform large scale studies in inexpensive modern computer platforms such as a small Pentium III cluster. The densityfunctional theory is the most accm-ate method [27, 41-42], and unlike the above mentioned methods, less insight, or educated guess, is required in the calculations. The density-functional calculations can also provide information for the less demanding tight-binding calculations, and is useful for vaUdating the results obtained from other studies. However, it demands much greater computer resources than tight-binding methods. Recently, large scale densityfimctional studies, such as the complex magnetic properties [43], are performed with the massively parallel processing computers. In the following sections, we first discuss the multi-band tight-binding Hamiltonian and the impiurity averaged Green's function. Then we will explain the calculation of the conductivity using the Kubo formula and the vertex correction needed for current conservation. In the last section, we discuss how to obtain PAGE 28 20 tight-binding and impurity parameters corresponding to experimental situation and compare some of the theoretical results with experiments. 2.1 Modeling of the Band Structure In this section, we discuss the modeling of the band structure using a tightbinding Hamiltonian. Then we explain how to model impiurities in samples and how to obtain the impurity averaged Green's function, which contains the band structure information. 2.1.1 Multi-band Tight-binding Hamiltonian . The general form of a tight-binding Hamiltonian is given by H= H{r,r';j,j';a,a')clj,Cr'ry, .. (2.1) where clj^ and Crj> are creation and annihilation operators of an orbital j with spin a at site r. Depending on the crystal structure and the details required in the model, the hoping terms, H{r,T'; jj';a,a'), are set to be zero when r and r' are third neighbors or fiu-ther. However, in bcc crystals, the third neighbor is essential. In modeling the transition metals, it is often important to keep all the orbitals in the outermost shell: s,p,d. The coupling between the majority and minority spins H{r,r';j,j';^,i) contributes to spin-flip scattering. For example, to include the spin-orbit coupling, a term Â• S is added to the Hamiltonian [44]. For transports in transition metals samples we consider, the spin-diffusion length is much larger than the layer thick- PAGE 29 21 ness [28-29], therefore spin-flip is neglected in most calculations. In this case, the two spin species axe assumed to be independent and with decoupled Hamiltonians, i/(t, T) ^nd H{l,i), and can be solved separately. Unless otherwise specified, the formulation is this work can be applied to the full Hamiltonian as well as the Hamiltonian for individual spins. 2.1.2 Modeling Impurities Scattering sources exist in all samples including in experiments of heterostructiu-es. If we do not include these effects in our calculations, the resistivity computed would be at least one or two magnitudes less than those observed in experiments. There are three major source of scattering: disorder in the bulk of the sample, impurity atoms usually due to interface diffusion, and phonon scattering. Figure 2.1 shows a clean multilayer lattice and the same lattice with interface diffusion, where 10% of each monolayer diffuses into each of the two monolayers above and below. One can think of an impiurity atom in a host as a new heterostructure. However, as we will explain in the next subsection, the advantage of treating it as an impurity is that simplification is often possible. For example, in an Co/Cu interface, we would treat the Co atoms diffused into the Cu layers as Co impurity in Cu host, and the Cu atoms diffused into the Co layers as Cu impurity in Co host. The impurity potential of an impurity B in an host A is defined as AÂ«(r,r') = HB{r,v')-HA{r,r') (2.2) PAGE 30 22 where Hb{t, r') and Ha{t, r') are the hoping between r and r' for the perturbed and unperturbed system. The most important contribution to Aw is usually the on-site term. As a perturbation, we found it sufficient to model this term as diagonal in all indices. In other words, Au{r,T';j,j';a,a') = Au{r,T-j,j;a,a)5rr'Sjj'S^a'. (2.3) Figure 2.1: Illustration of a multilayer without and with interface diffusion. Without the interface diffusion (a), the multilayer has translational symmetry in the in-plane directions. In (b), 10% of each monolayer diffuse into each of the two monolayers above and below. By adding interface diffusion, the translational symmetry is broken. Two types of disorders are included in our model: the substitutional disorder which is often caused by interface diffusion and the structural disorder which are randomly located in the sample. To model the substitutional disorder in the 3d transition metals, we set the substitutional disorder parameter PAGE 31 23 Awsub(r, r;j,j; a,a) = 0 for j Â— s,Px,Py,Pz because the most important difference between the tight-binding parameters in different elements is the on-site potential of the d-orbitals. The difference among the d-orbital parameters are small; therefore, it is often a good approximation to set Attsub(r> r; ifi; cr) = AMd(r, cr) for j = dxy,dyz,dzx,di,d2, such that two parameters, AMd(r, j) and Awd(r, J,), parameterize the impurity at site r. To model structural disorders, we use impurity potential of the form Au(r, r'; j, f; a, a') = AustSrr'SjfS^a' , (2.4) where Aitgt is a single parameter that characterizes the bulk of the sample. In this work, we do not model the phonon scattering. However, we can increase the structural disorder parameter for the bulk to get an idea of the effect of increasing temperature. 2.1.3 Impurity Averaged Green's Functions Once the Hamiltonian and the impurity parameters are given, we are ready to obtain the band structure information by either solving the eigenvalue problem or the Green's fvmction. Here we explain the impurity average Green's function approach which is especially useful in treating systems with impurity scattering. The impurity average Green's function is the ensemble average of the Green's function for systems with different microscopic impmrity configiurations in the same macroscopic impurity distribution. In many systems, geometric symmetries are broken by the impiu-ity atoms. For example, in the multilayer structiire in Figmre PAGE 32 24 2.1(b), the inplane translational symmetries in the unperturbed system (Figure 2.1(a)) are broken by the impurities due to interlayer diffusion. However, the impurity probabihty distributions in each site within the same layer are the same. Therefore, after averaging over the impurity ensemble, the system recovers the translational symmetry. Also, the average of the Green's function is performed analytically. As a result, the impurity averaged Green's function approach can be used to study much larger systems than techniques that requires averaging over the impurity configuration numerically. We first consider the system unperturbed by impurities. The unperturbed, or bare, retarded/advanced Green's functions is given by {ujHÂ±ir])Go^'^ = 1, " . Â• (2.5) where the scalar u is the energy, t] is an infinitesimal small positive number, and 1 is the identity operator. The Green's function and the quantities obtained from the Green's function are a;-dependent, in the following, the u) label is suppressed for clarity. From the definition, the bare Green's functions satisfy : {G^^y = Gi^^. , (2.6) Since the tight-binding parameters we use are real, the Hamiltonian is real^ in the real space representation. In other words, H'^ = W Â— H, which leads to This is not true in the Fourier space (2.7) PAGE 33 25 When impurity atoms are added to the system, one may try to compute the Green's function for system with impurities at their exact location. However, in most experiments, the exact location of impurity atoms are unknown, and even the probability distribution of the impurities are obtained by educated guesses. Moreover, the effect of the randomly distributed impurities are averaged out. Therefore, one has to compute the impurity averaged Green's fimction^, defined as {u-HSÂ«/^)G^/^ = 1. (2.8) where S^/"^ is the self-energy. For the on-site, diagonal impurity potential described in the last section, the self-energy is given by SÂ«/^(ri,r2) ^ ^(ri)Jr,r,G'Â«/^(ri,r2) (2.9) where the impurity parameter ^(ri) = {Au\r^)) (2.10) is the ensemble average of the square of the impurity potential. For convenience, we denote Equation 2.9 as Y^R/A ^ j;^^. (2.11) ^For convenience, we use G^/-^ instead of {G^/^) as the symbol for the impurity averaged Green's function. PAGE 34 26 In this formula, we have neglected the cross-diagrams, which corresponds to fourth order terms in the impurity perturbation [45-46]. Since the impurity average Green's function and the self-energy depend on each other, we have to solve Equations 2.8 and 2.9 self-consistently. Other types of scattering effects, such as those by the bulk disorder, can be modeled within the same framework. We can also write the Green's functions in Dyson's equation as G^M = G?/^ + gJ/'^E^/^G^/^. (2.12) The following properties follows from Equation 2.6: (GÂ«/^)t = G^/^ (EÂ«M)t = s^/fi. (2.13) In addition, in real space representation, , . (SÂ«/^)^ = SÂ«/^. '; (2.14) 2.2 Conductivity Calculation In this section, we explain how to compute the conductivity using the Kubo formula which is the theory of linear response. We first review different current related operators in the tight-binding model, which will be then be used in the PAGE 35 27 Kubo formula. We will also discuss the vertex correction needed to give a current conserving conductivity. 2.2.1 Formal Aspects of Current Operators In this and the next subsections, we discuss discuss four operators related to the current: the total current, the divergence of current density, the current through a siuface, and the current density operators. These operators will be useful for computing the conductivity and checking the ciurrent conservation. In this subsection, we will discuss these operators formally. In the next subsection, we express these operators in matrix form explicitly for an one dimensional Hamiltonian. Classically, the current density is defined as the charge density times the velocity. The total current, the divergence of current density, and the current through a surface can be obtained from integrating or differentiating the current density. In quantum mechanics, the current operator is defined using the Hamiltonian. However, in the tight-binding model, the current density is more difficult to define because of the discrete nature of the basis. The current density represented in a discrete basis, it is defined from the current from atomic sites r to r' divided by the cross-sectional area. In other words, there is only off diagonal element in the matrix representation. Unfike the classical current density, the cvurent density operator not a localized in space. On the other hand, the total current, the divergence of current density, and the current operators are easily found from the charge density operators though the current conservation equations in Table 2.1. We will use the current operator to define the current PAGE 36 28 name operator relation to Jr conservation equation total current J divergence Vr Jr ^a=l aa current / an ^ Â— ~'Â§i S(r-ro) n<0 ^rPr Table 2.1: Summary of current density related operators in the tight-binding model. The names, symbols, relations to the current density operator Jr, and conservation equations used are shown in the table. In the tight-binding model, we define the current operators using the conservation equations, and use them to define the current density operator. In the table, Vr is the volume per atom, pr is the charge density, P is the polarization operator, is the atomic spacing in the Qf-th direction. density operator, and use the operator for the total current and the divergence of the cmrrent density to check omdefinition. The relations between the operators discussed in this section are similar to those of the classical quantities; however, some details must be treated carefully. Table 2.1 summaries the current operators and how can they be obtained thru current conservation laws. The total cmrent operator can be obtained [47] from the polarization operator, which is defined as P = e ^ rc|:cr. (2.15) PAGE 37 29 The total current operator is given by J = ^ dt = i[H,P] = -ie^{T-T')H{r,r')clcr'. (2.16) r,r' Here, the Planck's constant divided by 27r, is set to equal to one. Since classically the total current is the volume integral of the current density, one may be tempted to define the current density operator Jr as proportional to X^r'(l/K')(r ~ r')[/f(r, r')4cr' Â— h.c.]. However, this definition violates current conservation because the divergence of this expression does not agree with the -dpr/dt. We find the operator Vr Â• Jr using the continuity equation K.J, = -^. (2.17) In tight-binding model, the density is given by Pr = ec\CrlVr, (2.18) where Vj. is the volume occupied by the atom at site r. The Schrodinger equation implies Vr-Jr = Â—ie Vr [Hclcr clcj-H] (2.19) PAGE 38 30 Â—le Y,[H{r\r)clcr-clcr'H{r,r')]. (2.20) It is useful to note that in matrix representation the only nonzero element of 4cr' is 1 at (r, r'). In other words, (4cr')(ri, = (^rri(^r'r2Therefore, we can write the element of the divergence of the current explicitly as (Vr Â• Jr) (ri,r2) i-ie/V,)[H{r,,r)Srr,-5r,rH{r,r2)] (2.21) The operator for a current through a plane containing point Tq and normal to a unit vector n can be found by (r-ro) n<0 = Â—te 5]F(ri,r2)4^Cr2 , 4cr 'ir2 (r-ro)-n<0 E E[^(ri>r)4iCr-/i.c.] (r-ro) n<0 ri -^^1 E E + E E ] {H{vur)cl^cr h.c] (r-ro)-n<0 (ri-ro)-n>0 (r-ro)-n<0 (ri-ro)-n<0/ = -i^ E E [H{r,,r)cl^cr-h.c.]. (2.22) (r-ro)-n<0 (ri-ro) n>0 We have used the fact that the second term in the fourth line is zero due to the symmetry between r and ri. Physically, the second term is zero because it contains only the hoping term that does not leave one side of the plane. The plane current operator obtained at the end contains only the hopping terms which bring current from one side of the plane to another. PAGE 39 31 Unlike the operators shown in the first column of Table 2.1, there is no simple equation that leads to the current density operator. In the continuous limit, the current density at r is given by the current though a infinitesimal siurface divided by the area. However, in the tight-binding model, we have to discretize the length scale to the atomic scale. As an example, in one dimension the current density operator in second quantized form is given by = -ia/K [H{xi,X2)cl^c^^h.c.]. " (2.23) Il>I>I2 One can check that the divergence, defined in one dimension as {Jx+i Â— Jx)/0', agrees with Equation 2.27, and J^Vx equals the total cm-rent J. This definition can be extended to higher dimensional lattices. For example, in a rectangular lattice with only nearest neighbors, the two component of the current density is given by = -iai/Vr [H{ri,T2)cl^Cr^h.c], xi>x>x2,yi=y2=y Jr' = -ia^/Vr Yl [H{rur2)clcr,-h.c.], (2.24) yi>y>y2,xi=x2=x where ai and 02 are the lattice spacings in the xand ^-directions respectively. 2.2.2 Current Operators in One Dimension In this subsection, we consider an one dimensional Hamiltonian and express the corresponding current operators in matrix form exphcitly. We can better un- PAGE 40 32 derstand the properties of the current operators by studying this one dimensional example. Consider a one dimensional Hamiltonian, H, with up to second neighbor hoping terms^. In matrix form, the Hamiltonian can be written as The Hamiltonian in matrix from is given by / H = \ H{b,3) H{5A) H{5,5) H{5,6) H{5J) if (6, 4) H{6,5) H{6,6) H{6,7) H{6,8) H{7,b) H{7,6) H{7,7) H{7,8) H{7,9) \ J (2.25) From Equation 2.16, the total current operator J can be found by Â—iea times the matrix \ 2H{5,3) H{5A) 0 -i/(5,6) 2//(6,4) H{6,5) 0 2i/(7,5) H{7,6) -2H{5,7) -H{6J) -2H{6,8) 0 -H{7,S) -2H{7,9) V where a is the lattice constant. ^The model with only first nearest neighbor is too simple to illustrate some of the properties. PAGE 41 33 The matrix form of the current density operator at, for example, x = 6, can be obtained from Equation 2.23 as Â—tea \ -i/(5,6) -F(5,7) i/(6,4) /f(6,5) ^(7,5) (2.26) where a is the lattice constant, and V is the volume occupied by a site, using this definition, one can see that the total current operator satisfies J = Ex JxV. We also note that in one dimension, the current operator is equals the current density times the area V/a. In other words, the current operator at x = 6 contain only the hoping terms which bring the current from one side of the plane x = 6" to the other side. By Equation 2.26, the divergence of the current density operator at x = 6 is given by dx x=6 a PAGE 42 34 -le T -H{6,4) -H{6,b) 0 H{7, 7) -H{6,7) -H{6,8) (2.27) This equation can be rewritten in the following form: dx -le x=6 V Hi4,6) 0 0 H{6,6) 0 0 H{7,e) H{8, 6) PAGE 43 35 / 0 0 //(6,4) Hi6,5) ff(6,6) H{6,7) H{6,8) 0 0 (2.28) which the same as Equation 2.20 in one dimension. If we sandwich this expression by two Green's functions or wavevectors, the first and the second terms of the expression can operate on the left and right respectively like the Hamiltonian. For example, the second matrix operated on the retarded Green's function on the right can be evaluated using Equation 2.8: This property of divergence of the current density is useful later in the disciission of current conservation. 5{x, x). (2.29) PAGE 44 36 2.2.3 Kubo Formula . . In the linear response theory, the non-local conductivity tensor is defined as the response function in the equation , ; ' Jr = 5:aÂ„'Er'K', (2.30) r' where Jr is the current density response at r, Er' is the applied electric field at r', and Vf' is the volume per atom. The impurity averaged conductivity tensor can be computed by the Kubo formula : 1 rr, 47r (2.31) In these equations, the ciurent density plus vertex correction T^^ is given by the Dyson's equations r;*^ = Jr + K^^^ (2.32) = 3r + vG^r^^G^ (2.33) = Jr + vG^JrG^ + vG^vG^3rG^G^ + .... (2.34) The other vertex corrections K^^, K^^, and K^^ are obtained by replacing the corresponding indices in the equations above. The Kubo formula is represented in the Feynman diagram is Figure 2.2. Since, K^'^ is diagonal and F^^ is not, in PAGE 45 37 practice, it is more convenient to iterate on K^^ instead of F^^: (2.35) vG^'irG'^ + vG^'K^^G^ . (2.36) We also define r' Air y ^ ^ ^ \ (2.37) For a constant electric field, Jr = arE. Therefore Oj. is often compared to experiments. It can be shown that this conductivity is still correct even in non-uniform electric field because the result depends only on the potential drop across the sample. As seen in the bubble diagrams in Figure 2.2, the conductivity can be found by first dressing up the vertex in either current operator and then multiplying it by the other current operator, followed by taking the trace. Analogous to J, we define T^^ = ErTf^K and = Zr^^^Vr, which gives = 3 + K (2.38) = J + i;G^r^^G^ (2.39) J + ?;G^ JG^ + vG^ vG^JG^ G^ + (2.40) PAGE 46 38 Similarly, the iteration equation for K^^ is (2.41) (2.42) Therefore, we also have 2tRAqAj^qR _ Y^^G^J^G^ T^^G^J.G^ (2.43) As we will discuss in the next chapter, the two sets of equations. Equations 2.322.37 and 2.38-2.43, have different computational costs in different situations. In Chapter 3, we will explain which set of equations should be used in cvurentperpendicular-to-plane and current-in-plane calculations. -IFigure 2.2: Feynman diagrams for the average conductivity, (a) The open bubble, which corresponds to the conductivity without the vertex correction, (b) The conductivity contribution from the vertex correction. PAGE 47 39 2.2.4 Current Conservation For a system without impiirity average, the conductivity is given by the the Kubo formula without the vertex correction, or the open bubble (Figure 2.2(a)). However, when averaged over the impurity configurations, the open bubble does not satisfies current conservation. This is because using the impurity averaged Green's function to compute the conductivity in the open bubble, one neglects the correlation terms such as those in Figure 2.2(b). These terms are required to give a conserving current. In the following, we show the terms in the vertex correction shown in Equation 2.34 will give a conserving current. In other words, we will show that the conductivity obtained from Equation 2.37 satisfies Vr-Jr = Vr-(7rE-0. (2.44) Since the electric field is arbitrary, we need to show Vr-(Tr = 0. (2.45) Since the current Jr' does not depend on r, we can first find the divergence of the vertex G^FrG^, and then multiply the divergence with the Jr', followed by taking the trace. The divergence of the zeroth order term of the AR vertex in Equation 2.34 is Vr-(G^JrG'Â«) = G^{Vr-Jr)G''. (2.46) PAGE 48 In the one dimensional case, this is simply a ^ ' a 40 (2.47) By current conservation (Equation 2.19), we have (Vr Â• Jr) = {-ie/Vr)[{G^H) 4crG^ G^lcr [HG"") = (-ie/Vr)[(-l + G^uj G^S'^)4crG^ G\lcr{-1 + uG"" S^G^)] = {-ielVr)[G^clcr 4crG^ G-^E^cJcrG^ + G^4crE^G^]. (2.48) Here, we have used the Schrodinger equation (2.8) to evaluate HG^/^ as (2.49) For the divergence of the first order term of the AR vertex (vertex with one rung) in Equation 2.34, we have Vr-(G^i;(G^JrG^)G^) G^'t;(G^(Vr Â• Jr)G^)G^ G^ vG^clcr GÂ« G^ vclcrG'' G^ -G^ vG^i:\lcrG'' G^ + G^ vG^lcr^'^G'' G^ (2.50) '5 PAGE 49 41 The first term of the above equation can be rewritten using the following equation: (t;G^4cr)(ri,r2) = v{ri)Â§r,,v2{G^4cr.){ri,r2) = 'y(ri)5ri,r2G^(ri,r)(^r,r2 = S^(ri,r)5r,r2 = (S^4cr)(ri,r2). (2.51) We can see the first term of Equation 2.50 cancels the third term of Equation 2.48. Similarly, the second term of Equation 2.50 cancels the last term of Equation 2.48. The same argument can be applied for higher order terms to show the last two terms of the divergence of the (n + l)-th order vertex cancel the last two terms of the divergence of the n-th order vertex. As a result, the only terms left of the divergence of the vertex are the first two terms of Equation 2.48. The divergence of the conductivity in Equation 2.37 is proportional to the sum of the following terms 2IV(JG^Vr Â• (r;^^)G^) = {-2ie/Vr)[Kf 4], -Tr(JGÂ«Vr Â• (rr^Â«)G^) = {ze/V,)K-4], -Tv{3G^Vr-{T^^)G^) = {ie/Vr)[Kf 4], (2.52) where Kf/^ = TV(JctcrG^/^) = Tr(4crG^/^J), kJ/^ ^ Tr(JGÂ«/^4cr), PAGE 50 (2.53) In the real space representation, = H, which impUes (4crGÂ«/^J)^ = -JGÂ«/^4cr. (2.54) Since transpose of a matrix does not affect its trace, we have /.f/^ (2.55) Unhke the continuous case [46], the divergences of the individual dressed bubbles in the tight-binding model are non-zero. However, they cancel each other so that the divergence of Gr obtained from the Kubo formula is zero. Therefore the vertex corrected conductivity satisfies current conservation. In this work, we only use on-site scattering potentials, the self-energy E^/'^ is also on-site. Therefore, E^/^4cr = 4crS^/^. (2.56) This means that in Equation 2.48 the last two term for divergence of the AA and RR vertex cancels each other out, and no vertex correction is needed for the AA and RR vertex. The Kubo formula for on-site scattering potentials can be simplified as ar = ^Tr[2JG'^(Jr + K,^^)G'^-JGÂ«JrG'^JG^J.G^], (2.57) PAGE 51 43 When the conductivity is calculated in nanostructures, for example in the current-perpendicular-to-plane (CPP) magnetic multi-layers, we find the contribution from the vertex corrections can be as large as 20% of the total conductivity. Therefore, the vertex correction should not be neglected. 2.3 Realistic Parameters In the earlier part of this chapter, we explained how to obtain the conductivity once the Hamiltonian and the impurity potential are given. In this section, we explain how to obtain the tight-binding parameters and the impurity potential, which are essential to modeling the electron transport accurately. We will also compare some theoretical results with experiments. 2.3.1 Tight-Binding Parameters The tight-binding parameters of an element is obtained from fits to the bulk band structure of density-functional calculations [48]. To do this, we compute the eigenvalues for a set of k-points in the irreducible zone with the density-functional theory, and also with the tight-binding model. The tight-binding parameter are then adjusted to minimize the error EENf(k)-^f^(k)P, (2.58) k j where ujj^{k), and u;?*^'^(k) are respectively the tight-binding theory eigenvalue and density-functional theory eigenvalue of the j-th orbital of a kvector. As shown PAGE 52 44 by Papaconstopoulos [48], tight-binding parameters obtained in this way give the correct band structure for bulk systems. For heterostructures, we do not obtain the tight-binding parameters by fitting to the band structure of the heterostructure because the density-functional theory for complicated systems can be expensive computationally. Since our goal is to model the system efficiently but still give reliable predictions, we use a less costly approach. The band structures for 3d transition metals are very similar, except for the offset of the c?-bands, which are different in each metal. To model 3c? transition metal heterostructures, we use the bulk parameters in the bulk and the interface. The on-site energy of each metal is shifted so that the Fermi-level level is zero. In this case, the Fermi level for different metal aligns. The hoping parameters between the atoms of two difference element are set as the average of their bulk values. 2.3,2 Impurity Potential for Substitutional Disorder As explained in Section 2.1.2, the substitutional disorder is caused by impm-ity atoms usually from interface diffusion. The impmrity potential for the atom can be obtained from density-functional fits. To the first order, one may use the difference in the bulk on-site parameters between the impurity and the host as the substitutional disorder potential, AusuhHowever, since in transition metals the on-site energy of the impurity B are often affected the host A, it is more accurate to find the impurity potential by fitting the tight-binding model to the density functional theory for a small unit cell such as A3B1. Again, the error between the tight-binding and density-functional eigenvalues given in Equation PAGE 53 45 2.58 is minimized. Since in the transition metals the most important difference between different elements is the d-orbital on-site energy, all the parameters except for the d-orbital on-site energy of the impurity can be set as the parameter of the host obtained from the bulk fit. As a result, the fitting for A3B1 is just an oneparameter fit for each spin, where the only varying parameter is the d-orbital on-site energy. The validity of the parameters are then confirmed in small unit cells by comparing the density of states computed by the tight-binding theory and the densityfunctional theory. For example, parameters for Fe as impurities in fee Ni host are obtained from fits to fee NisFe band structures. As a check, parameters obtained from the fit are used to calculate the tight-binding band structure of fee NisiFe, which agrees with the band structure obtained from the density-functional theory. This indicates that the model works well at least when impurity concentration is less than 25%. When the impurity concentration is large, or when the difference in hoping between the host and the impurity is large, the error is expected to increase. The other parameter in interface diffusion is the density of impurity atoms in the few atomic layers next to the impurity. Unfortunately this parameter depends the material and methods used to grow the sample in experiments. Only in a few transport experiments the sample is characterized well enough to give the details of the interfaces. Even in those case when the interfaces are well characterized, such as those by Cyrille et al. [6-7], only the long length scale (Â«10 nm) interface roughness is measured. Within our model, the interface roughness at this length scale should be considered in the geometry because the Fermi wavelength is an PAGE 54 46 order of magnitude smaller. On the other hand, the atomic scale interfaces roughness is related to the density of impurity atoms. In Chapter 4, we will explain how to set the density of the impurity atoms. 2.3.3 Impurity Potential for Structural Disorder The impurity potential for the bulk structural disorder, Awst, can be obtained by comparing the calculation of the resistivity and the resistivity for bulk samples. Usually, the samples in heterostructures are not as clean as a bulk sample. Therefore, we use the room temperature resistivity in bulk samples as a guide to setting the bulk disorder parameter. 2.3.4 Extension to Magnetic Alloys The method described in this section can also be extended to model transition magnet alloys, where the magnetic moment and hence the band structiure vary with the composition. In these materials, the localized d electrons are responsible for the magnetism. We can assume that the most important changes due to alloying are contained in the on-site energies, i/(r, r, d, d, t, t) and /f(r, r, d, d, |, |), of the spin-up (majority) and spin-down (minority) d electrons. In other words, the alloy is assumed to have site-independent hoping parameters and s and p on-site energies the same as the host metal, and site-dependent d on-site energies. This assumption will be justified by the agreement with local density calculations of supercells in different impurity concentrations. The major contributions to the on-site energies come from the atomic core potential plus the Coulomb energy due to the opposite spin. In the itinerant PAGE 55 ':. 47 electron model, it is more convenient to work with the total number of spin-up and spin-down d electrons, N}^ and A'^^, at each site r. Therefore, we rewrite the parameters in the form of F(r, r, rf, r, T) = U^ + U^Ni^d H{r,r,d,d,i,i) = U^ + U^NlÂ„ (2.59) where [7Â° and Up are the on-site parameter and effective Coulomb energy per pair at site r. Both A^J^ and A/^f^ in Equation 2.59 are obtained from the band structure, which in turn depends on .^^(r, r, d, d, j, t) and H{r, r, d, d, |, j). Thus, Equation 2.59 has to be solved self-consistently with the alloy band structure calculated using Hamiltonian of Equation 2.1. The alloy band structmre is computed with the coherent potential approximation (CPA) [49-50] is used, which is summarized in Appendix B. All parameters for the host atoms can be found from fits to the local density calculations of the pure metal [48]. The only two parameters left, C/j^p and 11^^^ of impm-ity atoms, are obtained from fits to the local density calculations of supercells composed the two types of atoms. In Figm-e 2.3, we plotted the Slater-Pauling curve, which shows the magnetic moment against the average number of electrons, for 3d magnetic alloys. Both experimental data and tight-binding CPA calculations are shown. As we can see, the theory agrees not only with the main curve, but also the NiCr branch as well. The NiCr curve is very sensitive to the details of the band structure. In this simple model, only two parameters in addition to the parameter of the host are PAGE 56 48 included; therefore, there is some deviation in the NiCr curve. The agreement indicates that the tight-binding model described here is suitable for obtaining the magnetic properties and band structure properties of transition magnetic metal heterostructures. CO E o i0) Q. C o E o E .0 o c O) CO E 7 8 9 10 number of valence electrons per atom 11 Figm-e 2.3: The magnetic moment of alloys per atom versus the number of valence electrons per atom. The lines are the experimental data of FeV, NiFe, NiCu, and NiCr alloys as labeled. The symbols are from the CPA calculations of FeCr (triangles), NiCr (circles), NiCu (squares), and NiCr (crosses). The CPA calculation agrees well with the experiment in the main branch (NiFe and NiCu), and agrees reasonably well with the NiCr branch. The dip in the experiment curve in NiFe corresponds to the bcc-fcc phase transition. PAGE 57 CHAPTER 3 COMPUTATION PROCEDURE In the last chapter, we discussed the formal aspects of modeUng magnetic nanostructnres and calculating the conductivity. We also discussed how to obtain parameters so that the models give reaUstic results that can be used to compare with experiments. However, if one writes computer codes directly from these equations, the resulting codes will not be fast enough for large scale studies. In particular, many of the equations involve 0{N^) operations such as the matrix multiplication of NM X NM dense matrix, where A*" is the niunber of atoms in a unit cell, and M is the number of orbitals per site. In order to speed up the calculation to tackle the tasks such as optimization and studying large systems, we have to reformulate the problem with special consideration about the computation speed. In this chapter, we discuss the computation procedure of the conductivity calculation. In the first section, we discuss the current-perpendicular-to-plane (CPP) conductivity computation, which demonstrates how to find the impurity averaged Green's function, solve for the vertex correction, and compute the conductivity. Although we use the CPP geometry as an example, the application to general geometry are similar. Most importantly, the tools that are developed here can be used in more general cases and the size-dependence of the computational cost are the same. As we will show, the size-dependence for computing the self-energy, the vertex correction, and the conductivity are 0{N), OiN"^), and 0{N), respectively. 49 PAGE 58 50 This is an 0{N) saving compared to total computation costs in the last Chapter. As a result, we can study more systems with a larger size. As we will see in Chapter 5, the speed of the code is crucial to optimization. In the second section, we will discuss the current-in-plane (CIP) geometry, which is also commonly studied in experiments. We will show that since the Hamiltonian has translational symmetry in the in-plane direction, the vertex correction is identical to zero. Therefore the CIP conductivity is easier to computer. However, the size-dependence of the overall computation remains 0{N^). In the third section, we briefly discuss the calculation of the conductivity tensor and the current-at-an-angle-to-plane (CAP) geometry, which will only require tools developed in the first two sections. In the last section of the chapter, we explain how we implement the algorithms discussed in this chapter. We shall use some new C-|Â— Imatrix hbraries, which will call the FORTRAN libraries BLAS (Basic Linear Algebra Subprograms) [51] and LAPACK [52] for the actual computation. We will also use an in-house blockmatrix library to implement most algorithms, and explain the advantage of using this new tools for coding. 3.1 Current-Perpendicular-to-Plane Conductivity We discuss the computational aspects of the current-perpendicular-to-plane (CPP) conductivity calculation in this section. Let us consider a periodic multilayer system that has N mono-atomic layers in each period and can be generated by a unit cell of A'' atoms. We label the growth direction as z and the in-plane directions as x and y. Since in the imit cell each atom belongs to a different mono- PAGE 59 ' . 51 layer, we label the an atom in the unit cell as z = 1,2, . . . , N. An example of this is a multilayer grown in the fcc(lll) direction. Each elements of a matrix is an m x m block matrix, where m is the nimiber of orbitals per site. The current is in the ^-direction. The equations in this section are intended for direct implementation, therefore, the computational cost of major equations are given. 3.1.1 Impurity Averaged Green's Function Using the model discussed in Chapter 2, the self-energy is given by : SÂ«(ri,r2) t;(ri)C;^(ri,ri)(J(ri,r2) (3.1) We only need to calculated the real space function for sites within the unit cell. Since S is diagonal, it can be stored in a vector and the computation cost is N. In real space, the Green's function is given by G^={uj-HE^)-^ (3.2) However, since a periodic system is also infinite, we have to calculate the Green's function for the imit cell in the Fourier space: (3.3) PAGE 60 52 where k is a wavevector in the first Brillouin zone, i/k and are the Hamiltonian and self-energy in the k-space. Since the self-energy is on-site, we have = in the unit cell. The Hamiltonian in k-space can be found by i/k = Hr -^^(0, R)exp(Â— ik Â• R), where R is the position vector of neighbor unit cells. In the following, we consider up to second neighbor hoping in real space. In the fcc(lll) lattice, up to the second nearest neighbors are included in the three monolayers described by ^ = 0, Â±1. Therefore, when we work in the k-space representation, the Hamiltonian is tridiagonal with periodic boundary condition. In Appendix A, we explain how the inverse can be found efficiently. Since ffk is tridiagonal with periodic boundary condition and E is diagonal, the cost for inverse is 0{N'^). However, since only diagonal elements of G are needed to find the self-energy, the cost can be cut to 0{N) in this step. The real space Green's function within the unit cell is then given by G^ = ^EGk (3.4) (Vk k The advanced self-energy and Green's function can be foimd by taking the hermitian conjugates of the results. The Equations 3.1 and 3.4 are solved self-consistently by iteration. If the system is finite instead, then we can work in the real space and the inverse algorithm for fixed boundary condition in Appendix A can be used; however, one has to consider the lead with special care. PAGE 61 53 3.1.2 Vertex Correction To find the conductivity, we have to compute not only the diagonal of the Green's function, but also the full Green's function. In Appendix A, we explain how to compute the Green's function in 0{N'^) for both periodic and fixed boundary conditions. In the k-space, the current density operator at site z in the unit cell is given byi Jzk{zi, Z2) = {-ia3/V,)Hi,(zi, Z2) a zi= z,Z2 = z + 1 ^ {ia3/Vz)Hk{zi, Z2) a zi = z + l,Z2-^ z = 0 otherwise, (3.5) where 03 is the monolayer distance. The total current is Jk = E'^^kK, (3.6) z which has 2A'' elements. As discussed in Subsection 2.2.3, both sets of equations. Equations 2.32-2.37 and 2.38-2.43, can be used to compute the conductivity. When we compute the vertex correction, both Equations 2.32 and 2.38 have the same computational costs. To see this, we notice from the definition of K^^ that it is non-local. Therefore, both K^^ and K^"^ have the same number of elements and are computed in the ^Here the coordinates are cyclic, for example, z = N +1 refers to the same site as 2 = 1. PAGE 62 54 same number of steps. Since is local, the cost for computing lix[GTGJz] in Equation 2.43 is 0{N), while the cost of Tr[GJGr,] in Equation 2.37 is 0{N'^). The vertex correction is given by the Fourier transform of the Dyson's equation (Equation 2.36): k 24,25 + E^(^i)^k (^1, z,)K^\z,, z,)G^{z,, z,) Z4 =^ ^i^l) E E[<^k (^1, 2:4)4(24, Z4 + 1)G^{Z4 + 1, 21) k 24 +G^izi,ZA + 1)4(24 + 1, z)G^iz4, 21)] +^(^i)E[G'k(^i' ^4)^^^(24,24)^1^(24, 21)] (3.7) 24 Again, K^^ = K^"^ in the unit cell because K^^ is diagonal. The first and second term cost 0{N'^), however, they are independent of the iteration. In the last term, for each 21, we need to do a sum over 24, thus, each iteration in this loop costs 0{N^). Since this step is inside the iteration loop, it is the most time-consuming part of the whole calculation. The cost for computing G^ to start the loop is also 0{N^). 3.1.3 Conductivity The conductivity is given by Fourier transforming Equation 2.43 as = JzET^ [24kGi^(4 + K^^)G^ 4kG^4Gfk J.uGUCt] -(3.8) PAGE 63 55 As explained in Chapter 2, is independent of 2, and hence we can choose any z the evaluate the conductivity. One should note that it only take 0{'N) to evaluate the traces. For example, the first two terms can be written as Tr +{G^J,y,G^){z',z' + l)Mz' + l,z'), ' ' (3.9) and V. J, Tr GUi^G^K^^] = Y.iGtJz^^G^){z\z')K^\z\z'l (3.10) where {GtJ.^.G^'\z'-^\,z') = {GtJ,^GZ'\z\z' + l) = {GUyG^^^)iz',z') = G^{z' + 1, z)J,y,{z, z + l)G^^^{z + 1, z') +G^{z' + 1,2 + 1) J,k(^ + 1, z)G^^^(z, z'), G^{z', z)J,^{z, z + 1)GZ'^{z + 1,2' + 1) +Gi^(2', 2 + 1) J,k(2 + 1, 2)Gf/^(2, 2' + 1), Gt{z\z)J,y,{z,z + l)G^/''{z + l,z') +G^{z',z + l)J,^{z + \,z)G^"'{z,z'). ' (3.11) Note that there is only one sum in each trace, therefore, the cost for each trace is 0{N). Further simplifications are possible by using the symmetry of the matrices. PAGE 64 for example, 56 Tr[(G^J,G^)(^', z' + l)J{z' + = Tr{[J{z', z' + J,GÂ«)(z' + 1, z')]^ = (Tr{(G^ J,GÂ«)(/ + 1, z')J{z\ z' + 1)})*. (3.12) Therefore, J,G^ J] = 2 ^ a^KG^ J,G'^)(z' + 1, z') z' + 1)}. (3.13) 3.1.4 Diagonal Impurity Potential Most of the discussions in this work allow the impurity potential to be offdiagonal in the orbital indices. However, for the diagonal impurity potential, v{z){m,n) = v{z){m)5{m,n), further simplification that reduces the cost of the second loop to negligible is are possible. In this case, since K^"^ is also diagonal in band index, we denote K^^{z){m) = K^^{z, z){m,m). In the following, we will show that the last term in Equation 3.7 can be written in the form of a matrixvector multiplication. Hence we can group the K^^ terms on both sides of Equation 3.7 and write this Dyson's equation in the form of il-q)K^^=p, (3.14) where g is a matrix and p is a vector. Equation 3.14 is a set of hnear equations and can be solved without iterations. To derive this equation, we notice that for PAGE 65 57 each zi, K'^{z,){m) = p{z,){m) and -^^^^ Y: E GÂ§{zÂ„Z2){m, n)K{z2){n)Gt{zÂ„ z,){n, m) k 22 = p{z,){m) + Y.q{z,,Z2){m,n)K''\z2){nl (3.15) 22 where ^^^^ E E[G^(^i> ^2) Jk(^2, ^2 + (^2 + 1, ^l) k 22 +Gk(-2l, 22 + 1)4(22 + 1, 22)Gk (22, 2l)](m, m) ^^^^^EE[Gk(^i>^2)Jk(22,22 + l)Gi^(22 + l,2i)](m,m) + h.C. k 22 E E Â»{[Gl^(2l, 22)4(22, 22 + 1)G^^(^2 + 1, 2i)](m, m)}. (3.16) g(zi,^2)(m,n) = ^ gÂ«(zi, Z2)(m, n)Gi^(^2, 2i)(n, m) ^ k y(2i)(m) El<^k(^i'^2)(m,n)|l (3.17) PAGE 66 58 Since the second term of p is the conjugate of the first term, the only term in p we need to compute is the first term. It can be calculated as follows : [G^{zu Z2)Mz2, Z2 + 1)G^{Z2 + 1, Zi)]{m, m) = E[<^k (2i> Z2)Mz2, ^2 + m')G^iz2 + 1, zi){m', m), m' (3.18) which does not involve any sum over positions. Since q is symmetric, we only need to calculate half of the element, and we can use the symmetric matrix solver for this problem. These two cost saving steps are important because they are 0{N^) calculations. When impiurity potential is not diagonal in the orbital indices, K^"^ and G do not commute and it is very costly to solve the linear equations, therefore K^"^ is solved by iteration. However, for the impurity potential used here, it is much faster to Equation 3.14 directly, which is a set of real and symmetric linear equations and can be solved with a cost of The cost for computing p and q are proportional to NkN'^M^ and NkN'^M'^, where Nk is the number of k-points. Therefore, most time is spent in evaluating p. 3.1.5 Total Computational Cost The computational cost is important to the large scale we discuss in the following Chapters, therefore it is given here. If A^^, are the number of iterations in the self-energy loop respectively, and Nk is the number of k-points, then the total cost for a conductivity calculation is roughly {4M^)Nk{6N^ + {l9Ns + 27) N) floating PAGE 67 point operations (flop), where the factor 4M^ is the cost for a complex 9x9 block matrix multipUcation. As an estimate, we take Nk ^ 1700, M = 9, A'' = 50, Ng Â— 10, the total cost is on the order of 10^^ flop. On a machine capable of about 10^ flop per second (1 GFlops), this should take a few minutes. In practice, we performed the calculations a computer cluster with ten 866MHz Pentium III CPUs, it takes about 5 to 10 minutes for a calculation when jV = 50. The algorithms presented in this chapter are CPU intensive. The memory requirement is 0{N^). For the computers used in this work, unless N is larger than about 500, one does not need to rearrange terms to optimize the memory usage. 3.2 Current-In-Plane Conductivity The current-in-plane (CIP) conductivity is a special case where the Hamiltonian has a translational symmetry in the direction of the current. Let us consider the same fcc(lll) lattice discussed in the last section and denote the direction of the current as x. The problem is greatly simplified because the vertex correction AT^^ is zero. This result can be shown by studying the first correction term of Equation 2.40: (^G^7G^)(x,x) = v{x)[J2G''{x,x + j)Jix+j,x + j + l)G^{x + j + l,x) j + ^ x + j + l)Jix + j + l,x+ j)G\x + 3, x). j (3.19) PAGE 68 60 In the real space representation, the Hamiltonian is real and has translational symmetry, hence and are symmetric, and J is antisymmetric. The second sum above can be rewritten as J2 G^i^ + i + 1, r)J{x + j,x + j + l)G^{x, X + j) J = G^^i^^ ^ i ^Ui^ -j-hx3)G\x i, x), (3.20) i where in the last line we have iised the translational symmetry separately for G^, J, and G^. By substitute / = j + 1, we can see the two sums in Equation 3.19 cancel exactly. As a result, all the correction terms to F^^ in Equation 2.40 are also zero. Since the vertex correction is zero, we can compute the conductivity corresponding to the open bubble by c^a; = Â— ^Tr[2JkG'j^Jj;kG^ JkG^JxkG^ JkG^JxkG^]. (3.21) Since J^k is a tridiagonal matrix in the unit cell, it has SN elements. Thus, the computational cost for each trace is 0(A''^). 3.3 Conductivity Tensor As explained in Chapter 2, the current density operator is not obvious to construct. It is especially difficult when the direction of cmrrent is at an angle to the crystal axes. Therefore, in computing the cmrent for systems such as the current- PAGE 69 at-an-angle-toplane (CAP) system, it is easier to compute the conductivity tensor and use the tensor equation J = crE. In a cubic bulk system, aÂ°'^ = Sa^cr^^However, in nanostructures, the symmetry is usually broken, therefore aÂ°'^ are in general different. The techniques for computing current with vertex correction is useful for computing the conductivity tensor. 3.4 Coding As shown in the above, it is natural to formulate our problem in block matrices and hence nattual to write our code using block matrix object classes. We have to implement the block matrix algorithms such as inverse and products discussed in previous sections to take advantage of the sparsity of the matrices. Therefore, it is not enough to just have a block matrix class and common algorithms. The block matrix class should help us develop new block matrix algorithm efficiently. Our goal is to translate mathematical equations involving block matrices into programming lines directly, while having performance close to Fortran codes. In this section, we discuss how this is achieved with C+-I-. First we use the Template Numerical Toolkit package to handle the operations of all dense matrices. To increase the speed, we link critical subroutines, for example, the matrix multiplication or inverse operation for dense blocks, to the Fortran/ Assembler coded native binaries BLAS [51] and LAPACK [52]. Next we develop a block matrix class that can use the TNT matrix class as elements. Using the new block matrix class, we develop codes to implement the algorithms described by previous sections. PAGE 70 62 In addition to these steps, our code is parallelized with PVM (Parallel Virtual Machine) at the level of k-points to use the all available computing resources of a heterogeneous network. The resulting code should scale well with the number of processors. 3.4.1 Template Numerical Toolkit (TNT) The template [53] feature in C++ language allows writing a single code for . : object classes or algorithms for different storage types. For example, the subroutine template PAGE 71 63 our platforms at the time we started coding.^ On the other hand, the TNT package could be compiled and linked with Fortran programs on our major platforms: Solaris and Linux. Consequently, we chose the Template Numerical Toolkit. The TNT object we use most are the Fortran_Matrix object, which stores a matrix using Fortran (column) convention, and the Vector object. They allow Fortran (one-based) indices such as A (2,3) . For example, the following lines define a 2 X 2 matrix A, vectors x and b. Fortran_Matrix< int > A(2,2, "12 14"); Vector< int > x(2, " 12"), b(2) ; b(l)=3; b(2)=4; The package also takes advantage of the function/operator overload feature of C-I-+ to make the matrix and vector operations intuitive. Basic operations like =, +, *, and Â« are overloaded. Continued from the above example, the next two lines computes and prints out the vector c = SA^x + b. Vector< int > c=3*A*A*x+b; cout Â« c Â« endl; Notice that the compiler links the three multiplications above to scalar-matrix multiplication, matrix-matrix multiplication, and matrixvector multiplication subroutines in order. Without operator overloading, the line that computes c would look like Vector< int > c(vecadd(vecmult(matmult(scalmult(3,A) ,A) ,x) ,b)) ; 2 According to the Blitz+-|website [57], Blitz++ can now be compiled by compilers like gcc and KAI C++. PAGE 72 64 3.4.2 Fortran performance with BLAS and LAPACK To achieve the performance close to that of a Fortran code, we overload the most often used and most numerically intense subroutines by linking them to BLAS and LAPACK. For example, we overload the matrix multiplications A*B by the BLAS matrix multiplication subroutines xGEMM, where x designates the data type of the matrices. This is done efficiently because the Fortraa_Matrix object stores matrix element using Fortran convention. To make our code easy to use and understand, we name commonly used functions similar to Matlab commands. For example, to find the inverse of a matrix, we use B=inv(A), which is linked to the LAPACK subroutines xGETRF (LUfactorization) followed by xGETRI (inverse of an LU-factored matrix). Since owe code spends most of the time doing matrix operations implemented with Fortran/ Assembler native binaries, the overall performance of our code is close to a pure Fortran code. On the other hand, the codes are easy to write and understand. 3.4.3 BlockMatrix class We design a BlockMatrix class to implement mathematical algorithms as directly as possible. Without such an object class, block matrix algorithms are difficult to code and debug. In the following, we explain which featines we want in a BlockMatrix class, and how these features enable us to implement comphcated block matrix equations without much effort. We want to design a BlockMatrix class as a matrix Fortran_Matrix objects. Using the Fortran_Matrix object, the elements (blocks) of our BlockMatrix class PAGE 73 65 will be intuitive and fast to operate. The storage should be efficient, which means that block matrices of different sparsity have different storage arrangements and different index functions. To do this, we define several subclasses such as the Full_BlockMatrix, PeriodicBand-BlockMatrix, and the Sparse_BlockMatrix. The index functions are defined within the subclasses to overwrite the virtual index functions defined in the parent class. The object should be intuitive to use. For example, one may obtain the current operator from a tight-binding Hamiltonian with the following code: for (int r=l; r<=N-l; r++) { J(r+l,r)=i*H(r+l,r); J(r,r+l)=dagger(J(r+l,r)) ; } Here, H(j , j+1) is a Fortran_Matrix itself, and also an element of the BlockMatrix H. In addition, we can access elements of the matrix H(j , by an expression like H(j,j+l)(ml,ni2). These featm-es allow us to implement block matrix algorithms efficiently and more importantly, debug in a shorter time. In the following, we show how the of Godfrin inverse method is implemented using the BlockMatrix class. We want to find G = where M is a tridiagonal block matrix. These matrices are stored as PeriodicBand-BlockMatrix, which is a subclass of BlockMatrix that stores a banded block matrix with periodic boundary conditions. For example, M is declared by PAGE 74 66 PeriodicBand_BlockMatrix< Fortran_Matrix< complex< double > > > M(N, KL, KD, KU, p, q) ; where N is the number of blocks, KL = KD = KU = 1 are the number of bands in the lower triangle, diagonal, and upper triangle or the matrix, and the block size is p X q. The Godfrin method shown in Equations A. 21 and A. 11 finds the inverse of a tridiagonal block matrix M in 0{N'^). To implement this algorithm, we first declare the subroutine, with input matrix M, and output matrix G and other variables used as temporary memory. Values such as the matrix size and block size are initialized, as shown below.^ template PAGE 75 The next code segment is directly translated from Equation A. 21. 67 X(N)=0.0; for (int j=N; j>=2; j~) { W=inv(X(j)-M(j,j)); C(j)=W*M(j,j-l); X(j-l)=neg(M(j-lJ)*C(j)); Y(1)=0.0; for (int j=l; j<=N-l; j++) { W=inv(Y(j)-M(j,j)); D(j)=W*M(j.j+l); Y(j+l)=neg(M(j+l.j)*D(j)); } Since in Equation A. 21, Wj ans Zj are not used later, the memory of W is reused. Notice that variables like M(j-l,j) and W are themselves matrices. Although M is tridiagonal, there is no need to worry about where M(j-l,j) is stored when writing this code. The details of the storage are handled in the index function of the BlockMatrix object class according to the sparsity of the M. What is left is straight forward translation of equations. This code segment should be transparent even to people without C++ experience, and it is therefore easy to debug. PAGE 76 The last part of the subroutine computes the Green's function using Equation A. 11. Again, this is a direct translation. for (int j=l; j<=N; j++) { G(j.j)=inv(M(j.j)-X(j)-Y(j)); for (int jl=j+l; jl<=N; G(jl, j)=C(jl)*G(jl-l, j) ; for (int jl=j-l; jl>=l; jl--) G(jl, j)=D(jl)*G(jl+l, j) ; } } The time spend on designing the BlockMatrix class should be less then the time we saved in coding and debugging other part of the code. Since the BlockMatrix object classes developed in this project are reusable, future implementations of mathematical algorithms involving block matrices should be as easy as the example shown above. PAGE 77 CHAPTER 4 APPLICATION TO GIANT MAGNETORESISTANCE In the previous chapter, we explained how to compute the transport properties in magnetic nanostructures. In this chapter, we apply the method to compute the current-perpendicular-to-piane giant magnetoresistance (CPP-GMR) [4-14] in magnetic multilayers. As explained before, a good theory has to agree with both the conductivity and the GMR. To date, we have found two sets of experiments on the Fe/Cr CPP-GMR with both published resistivity and GMR data, one is by Gijs et al. [4] and the other is by Cyrille et al. [6-7]. In the following, we first review these experiments. The we try to see if we can use theoretical models to reproduce the experimental data for reasonable values of the parameters. In order to understand the trends shown in the experiments, the magnitude of the computed resistivity be similar to the values found in the experiments. The changes in resistivity and GMR as a function of different parameters are much smaller in scale, they should be studied after the magnitude of the resistivity is understood. Using the model in Chapter 2, we are able to calculate both the resistivity and GMR values which axe close to Cyrille's data within reasonable parameters. However, the resistivity data of Gijs seem to be too small to fit in our theory unless very small disorder is used. We believe the required disorder is too small to be physical. 69 PAGE 78 70 4.1 Experiment of Gijs and Collaborators The Fe/Cr superlattice in the experiment Gijs et al. [4-5] are micro-structured pillars, shown schematically in Figure 1.2. The pillar in the middle in Figure 1.2(d) is the sample to be measm-ed. The samples sizes and measured conductivity and OMR are summarized in Table 1. Table 4.1: Summary of Gijs [Fe(30 A)/Cr(tcr)]ioo samples Cr thickness ^cr A Area GMR Psat(r = 0) 10 90 1.08 11.0 28 20 0.08 24.3 40 6 0.04 25.2 The data are presented in terms of GMR = Ro/Rsm 1, shown in Figure 4.1, and the quantity \J{RoRsa.t)RoA'^, shown in 4.1. Using these two quantities and the areas given in table 4.1, we can obtain the resistance Rq and i?sat, and hence po and psat. At low temperatiue, the resistivity of the samples is about 11 fiUcm for Cr layer thickness tcr = lOA, and is increased to about 25 nilcm for ^cr = 28A and 40 A. The GMR of the sample is 108% for ^cr = lOA and drops to 4-8% for tcr = 28A and 40 A. PAGE 79 71 As we shall see, the resistivity for both saturation field and zero field are very important to comparing with the theory. 120 I , T(K} Figure 4.1: Temperature dependence of CPP and CIP GMR of pillar structures by Gijs et al. [4]. 4.2 Experiment of Cyrille and Collaborators Cyrille et al. [6-7] studied the effect of the CPP-GMR as the number of bilayers N of the Fe/Cr superlattice increases. They found that the GMR increases as A'' increases. PAGE 80 72 0 L 1 1 1 I I 0 100 200 300 T(K) Figure 4.2: Temperatiire dependence ^J{Ro/?sat)-Ro^^ for in the data of Gijs et al. [4], where A is the area of the sample. This quantity is used in the ValetPert theory [16]. Figures 4.2 show the resistivity and the GMR of a series of [Fe(30 A)/Cr(12 A)]jv superlattices, where A'^ ranges from 10 to 60. In this figure, the resistivity for both saturation field, Psat, and zero field, po, is plotted as a function of the number of bilayers A'^. The magnitude of the GPP resistivity is about 35 fiQcm for for saturation field, and is about 40-50 fxi^cm for the zero field. The experiments also show that psat is roughly independent of A'' and po increases as A'' increases. This trend is not very clear on Figiu-e 4.2(a), probably due to sample variations. However, since po and psat are correlated, the trend is more PAGE 81 73 10 20 30 40 50 60 N Figure 4.3: CPP Resistivities and GMR in FeCr superlattice as a function of number of bilayers measured by Cyrille et al. [6-7]. In (a), the zero field resistivity po increases as the number of bilayers N increases while the satiuration field resistivity Psat remains roughly constant. In (b), the CPP-GMR increases from about 25% to 40% when N increases from about 17 to 40. This is a 60% change in the GMR. PAGE 82 74 clearly seen from the difference in resistivities which shows up in the GMR data in 4.2(b). If one concludes from the Figure 4.2(a) that Psat is roughly constant, then Po must be an increasing fimction of A^. By detail analysis on the property of the samples, the authors show that the long length scale (about 10 nm) roughness of the interfaces increases as A'^ increases, and suggest that the increase in po but not Psat is due to the increase in spin dependent scattering at the interface. Since we only include atomic scale interface roughness in our calculation, we would not give the same trend as the experiments. 4.3 Comparing Theory with Experiments In Figure 4.4, we have plotted both the data from Gjis et al. experiment and those from the Cyrille et al. experiments on the same graph. Cyrille's data appears to be on the y-axis because they are measured at low temperatures, which is required for the leads to be superconducting. As we can see, the Gjis' resistivity is about a factor of two larger than the corresponding data of Cyrille. Our calculations are at zero temperature so we do not explicitly have temperature effects; however, within our model, we have several disorder parameters. By increasing the bulk disorder parameter, we can increase scattering, which is one of the effects of raising the temperatmre. In Figure 4.3, we have plotted the the maximum and saturation resistivity for our calculations as the spin-independent bulk disorder varies. A moderate residual bulk disorder is also added to the system. The residual bulk disorder is chosen so that resistivity calculations of bulk Fe and Cr films agree with experiments. The results in Figiure 4.3 are close in magnitude to Cyrille's data and about twice larger than Gijs' data. PAGE 83 75 60l Â• 1 r 50% 100 200 300 400 500 600 T(K) Figure 4.4: Resistivity in Cyrille's [Fe(30 A)/Cr(12 A)]n samples and Gijs' [Fe(30 A)/Cr(10 A)]ioo sample plotted on the same graph. The saturation field (triangles) and zero field (squares) resistivity in Cyrille's experiment are close the the y-axis because the experiments is done at the Uquid helium temperature. The saturation field (pluses) and zero field (circles) resistivity in Gijs' sample are about 50% smaller than those in Cyrille's experiments. In order to fit Gijs' results, we have to calculate the resistivities using zero residual bulk disorder. The resistivities calculated this way are close to Gijs' data. However, the assmnption that the bulk of the sample is perfectly grown is unreasonable. There are other disorder parameters related to the interface that we can change. In this calculation, we assumed that the interface thickness to be 4 monolayers, and the interface mixing percentage to be 20%. The resistivities can be reduced PAGE 84 76 600 spin-independent bulk disorder x 4 x 10 Figure 4.5: Calculated resistivity in [Fe(28.5 A)/Cr(11.5 A)]oo for moderate residual bulk disorder. Both the resistivity for parallel (squares) and anti-parallel (crosses) increase with the bulk disorder, however, the parallel increases faster, therefore the GMR will decrease with the bulk disorder. The theory data is close in magnitude to Cyrille's data and about twice larger than Gijs' data. by reducing the interface disorder. However, we doubt that the interface can be much sharper than 4 monolayers and have less interface mixing. A better way to look at the differences among the data is to plot the change in resistivity against the saturation resistivity, as shown in Figiure 4.3. For a fixed temperature and sample, there are two independent data, the resistivity at zero field and saturation field. By plotting functions of these two independent variables against each other, we can identify trends and remove the uncertainty in temperature dependence in the calculation. As seen in Figure 4.3, Gijs' data and PAGE 85 77 the theory have the same trend, due to a change in temperature or a change in bulk disorder, Cyrille's data represents the change in number of bilayers and has a different trend. 25 20 E o CO V) Q. <10 Â» o 8 Gijs Theory Schuller 0. 0 10 20 Psat 40 50 60 Figure 4.6: Change in resistivity against saturation resistivity. The data of Gijs (crosses) and the theory (squares) have the same trend; while the data of Cyrille (circles) has a different trend because it is done by varying the number of bilayers. Another advantage of plotting the data this way is when we compare two curves on the graph, we are comparing parameters other than the bulk disorders, which are the most difficult to model because they are usually not directly measured in the experiments. With the most uncertain parameter eliminated, a good theory should fit the experimental ciurve well. PAGE 86 78 As seen in figures 4.4 and 4.3, Gijs' data has a lower resistivity than both the theory and the experiment. The trends in Figure 4.3 do look similar, but there is a large offset between the curves. Fred Sharifi suggested one possibility is that there may be other conduction paths that allow electrons to by-pass the Fe/Cr superlattice. In Figure 4.3 we can obtain Gijs' curve by adding a shunt resistor to the calculated superlattice model. Superlattice Figure 4.7: Parallel resistor model. The effective resistivity is found by 1/peff = 1/Piatt + 1/Pshunt5 where Piatt and Pshunt are the resistivity of the superlattice and the shunt resistor, respectively. In Figure 4.3, we plotted the the change in resistivity, po Â— Psat) against the saturation resistivity for both Gijs' data and our calculation with and without the shunt resistor. We found that the calculation with a temperature-independent, 80 ju!^cm shunt resistor agrees well with Gijs' result. Therefore, Gijs data can be explained if in the experiment, some electrons are allowed to by-pass the superlattice and thereby reducing the measured resistivity. PAGE 87 79 18 16 14 E12 o Sio (5 in Qo .^6 2I Â« Gijs o Theory : shunt Theory : no shunt 10 20 30 Psat (^^" C"^) 40 50 Figure 4.8: Change in resistivity, po Â— psat, against the saturation resistivity, PsatThe theory data without the shunt resistor (circles) are different from that of Gijs [Fe(30 A)/Cr(10 A)] (crosses), however, the theory data with the shunt resistor (squares) agrees with Gijs results very well. In order to fit the data, the resistivity of the additional path has to be almost independent of temperature and non-magnetic. It is not clear what caused the extra path. One possibility is that pin-holes are formed in the insulating polyimide, which allow conduction paths of Au to be deposited parallel to the sample. After the polyimide are washed away, the Au conduction paths are left on the sample. Another possibility is that the conduction paths may come from impurities at the edge of the superlattice. In Figure 4.9, we compare the theory with Gijs' [Fe(30A)/Cr(28A)] data. Again, with the shimt resistor, the theory agrees with the experiment very well. PAGE 88 80 The shunt used in this case has a resistivity of 75 /il2cm, this value is very close to the one obtained for the Fe(30 A)/Cr(10 A) sample. This suggests that the two paths are made of similar materials. By using our superlattice model with an extra conduction path way, we can explain both resistivity and GMR data for two different samples in Gijs' experiment. 18 16 14 Gijs Theory : shunt Theory : no shunt E12 o CO in Cl .^6 420 0 10 20 , ^30 , Psat 40 50 Figure 4.9: Change in resistivity, po Psat, against the saturation resistivity, psatThe theory data without the shunt resistor (circles) are different from that of Gijs [Fe(30 A)/Cr(28 A)] data, (crosses), however, the theory data with the shunt resistor (squares) agrees with Gijs results very well. PAGE 89 81 4.4 Other Effects on the GMR In this Chapter, we have assumed the adjacent magnetic layers to be completely anti-parallel at zero magnetic field. In experiments, the adjacent magnetic layers may not be completely anti-parallel due to the interlayer coupling [58-62]. The interlayer coupling controls the angle between the orientation of magnetization of adjacent layers at zero field. In Figure 4.10, the experiment of Parkin shows the saturation field oscillates as the thickness of the non-magnetic spacer layer varies. This effect also means that the magnetic coupling oscillates as a function of the layer thickness. The origin of the oscillation is the RKKY [58-59, 63-64] interaction. The GMR is also shown to oscillate with the thickness in the experiment of the Michigan State University group [9-10]. The effect is due to the GMR affected by how the magnetic layers align. Spacer Layer Thickness (A) Figure 4.10: Dependence of saturation field on spacer-layer thickness for families of Co/V, Co/Mo, and Co/Rh multilayers [59]. PAGE 90 82 200 1S0 s e S 100 o Â« c o Â« S FaiSnm) lCo< 1.5nm) / Cu(cl)l 0 CPP(Ho) Â« CPP (Hm) o CIP (Hm) 5.0 7.5 ^Cu Figure 4.11: The magnetoresistance of Co/Cu layers as a function of the thickness of the space layer (Cu) for current perpendicular to the plane of layers (CPP) [9-10]. Another complication is due to the fact that the magnetic layers are adjacent to non-magnetic layers or the vacuum. In such case, the magnetization in both the magnetic layer and the non-magnetic layer can be affected. More detailed information can be reviewed by studying the spin-density, such as those by Ohnishi et al. [42]. For example, Figm-e 4.12 shows that the spin density on the Fe layers extends to the vacuum, and there are negative spin density regions. However, since the computational cost in this effect is high, the GMR is usually computed PAGE 91 83 by assuming the magnetic layers are perfectly aligned in the saturation field and perfectly anti-aligned in the zero field. SPIN DENSITY Fe (OOl) Figure 4.12: Self-consistent spin-density map of seven-layer Fe(OOl) in unit of 0.0001 a.u. on the (110) plane. Each contour line differs by a factor of 2. Dashed lines indicate negative spin density [42]. PAGE 92 CHAPTER 5 OPTIMIZATION OF GMR STRUCTURES In this chapter, we discuss the search for the optimal atomic configuration for the giant magnetoresistance. In previous chapters, we have discussed how to model the transport properties of magnetic nanostructures and compared the theory to experiments. Given a good model, we should be able to make predictions as well. In particular, what atomic configuration gives the highest GMR? This question is difficult because the configuration space is large. Take magnetic multi-layers made from Sd transition metals as an example, suppose we were to make a multilayer that consist 50 monolayers, and we can choose to use among Co, Ni, or Cu atoms in each monolayer, there are roughly 7 x 10^^ possible configurations. It is hard to tell which configuration gives the highest GMR. Although theories are capable of computing the GMR for any of the configuration, it is impossible to exhaust every possibility. The conventional way of solving such a problem relies on physical insights, innovative ideas, and often trial-and-error combined with good luck. For a complicated function in the large configm-ation space, such as the GMR, explore within the subspace is limited to one that is simple enough to be understood with physical insights. In recent years, the combinatorics [65] method has been a useful tool in searching for new systems for interesting physical properties. For example, by systematically studying samples on a two-dimensional grid of gradual gradient of compositions, one can search though the configuration space at a much faster pace than was previously possible. However, trial-and-error 84 PAGE 93 85 methods are unlikely to locate solutions in, for example, the 50-dimensional space mentioned. In order to search for the solution in a high dimensional space, the computer program must not search the configuration space randomly, but also make use of the local information. Optimization methods such as simulated annealing and genetic algorithms can be used. For example, using simulated annealing, Franceschetti and Zunger [37] are able to find the largest bandgap configurations in Alo.25Gao.75As alloys in a unit cell with 128 atoms by searching for about 10^ steps, which is much smaller than the size of the configuration space. As discussed in Chapters 2 and 3, the GMR is complicated function of the atomic configuration and the disorder parameters. In the following chapter, we first review the simulated annealing method. Then we will apply the method to find out which configuration has the highest GMR, and discuss the finding and its implication. It is hoped that experimentally one can grow samples according to the optimal GMR configuration with non-equilibrium growth techniques such as the molecular beam epitaxy (MBE). 5.1 Inverse Problems Let / be a function of the atomic configuration C. We may want to find the minimum (or maximum) of the function in the configuration space /rain = niin/(C), (5.1) PAGE 94 86 and the corresponding configuration CminIn other situations, we may want to find the configuration that satisfies the target condition /(C) = /*^set^ which can be found be minimizing of the absolute value of difference. More complicated requirements, such as prescribing multiple physical property, can also be formulated as the minimization of the function [37] 0{C) = E^.l/i(^)-/r''l' (5-2) j where Wj is the weight assigned to the distance of property fj and the target /j^"^^^*. This problem belongs to a large class of "inverse problems" . These problems axe difficult because the functions in a high dimensional configuration space are complicated, and intuitive pictures are often not available. Finding the highest GMR configuration is an example of an inverse problems. In order to solve the inverse problem numerically, one must evaluate the function, or solve the "direct problem," many times. In practice, the direct problem itself is often difficult to solve. For example, it costs about 5 x 10^^ floating point operations to find the GMR of a superlattice with a unit cell of 50 atoms. To solve inverse problems, numerical optimization algorithms such as simulated annealing [32] and genetic algorithms [33-34] are used. These optimization techniques have been used to solve physics problems such as finding the equilibrium structure by minimizing the energy. However, the optimization describe in this chapter is rather different in nature. In solving the inverse problem, the solution obtained satisfies a set of conditions we prescribe, not the condition set by natiire. This allow us to tailor-make a configuration with the desired properties rather PAGE 95 than settle for what is naturally available. The optimized configuration may be realized by growth techniques that can build structures monolayer by monolayer, or even atom by atom. In addition, the configiiration space is very different from those used in the finding equilibrium structmes. Here, we can change the chemical composition in the configuration, which is not a thermodynamical process. 5.2 Simulated Annealing In this section, we explain how to maximize a function by the simulated annealing algorithm. The flow chart of the algorithm is shown in Figiu-e 5.1. To maximize a function / as a function of the configuration C, an initial configuration C is randomly generated and an initial "annealing temperatiue" T is set. The value of f{C) is then evaluated. A new configuration Ctriai is then generated by a Monte Carlo move from the previous configmation, and the value of /(Ctriai) is also evaluated. If the functional value of / for the new configuration is higher than that for the previous configuration, it is accepted. On the other hand, if the new GMR is lower, it is only accepted with a probabiUty of exp( (/(Ctriai) /(C))/T). The process is continued as the annealing temperature is gradually lowered until there is no further change in the configuration and the final configuration is a near-global maximum. This process resembles the annealing process in statistical mechanics, where the total energy is minimized when a system is cooled slowly. In simulated annealing, the value of / is used to decide whether to keep or discard the trial configuration. Compared with random walks routines, simulated annealing causes the system to spend more steps near the configurations possessing large values of /. Simulated armeahng is easy to implement, and it is typical that PAGE 96 88 the number of configuration tested in simulated annealing be much smaller than the size of the configiuration space. 6 Set initial configurati and pseudo temperature ioncX ture tJ Monte Carlo move C to Ctna! Find Af=-J(Ctna!)-f(C) Accept Gnal 5 Accept Ctnai with a probability P = exp( Af/T ) aowly decrease r to 0 e.g. T=Toexp(-tv'no) Yes >^ End ^ Figure 5.1: Flow Chart for Simulated Annealing. The flow chart for using simulated annealing to search for the optimal configm-ation C, that maximize a function /, is shown. 5.3 Accuracy and Speed Requirements on Model Searching for the optimal structure using optimization techniques requires a delicate balance between acciuacy and speed. On one hand, the model has to be accurate enough so that the optimal solution is reliable. The model has to PAGE 97 89 capture the physics and to accommodate a wide range of configiirations. The level of detail in the model varies from problem to problem. For optimization of the transport properties for nanostructiires, the model has to include important effects of the band structure of the materials, the geometry of the nanostructures, and the scattering effects caused by impiirity atoms and structural disorders. All of these factors significantly affect electron transport. This requirement eliminates methods containing less microscopic detail than the multi-band tight-binding model that we discussed in previous chapters. On the other hand, calculations of the transport properties in this model have to be fast enough so that it can finish in reasonable time. The number of annealing steps is expected to be of the order of 10^ to 10^. To finish the optimization within a month, one has to be able to compute the function for a configuration in about 4 to 40 minutes. This is very demanding even for running tight-binding codes on a computer cluster. At first glance, the matrix operations of finding the Green's function, vertex corrections, and conductivity in Chapter 2 are are all 0{N^) operations, where A'^ is the number of atoms in the unit cell. In addition, the self-energy and the vertex corrections are solved by self-consistent iterations. Performing optimization with a self-consistent 0{N^) calculation would take about 50 CPU-years, which is impractical. This is why we have developed the fast algorithm in Chapter 3. For A'' = 50, this algorithm is more than 100 times faster than the 0{N^) algorithm. Running on a Linux cluster with ten 866MHz Pentium III CPUs, it takes about 20 minutes for each OMR calculation and about two weeks for an optimization with 1000 annealing steps. PAGE 98 90 The computational cost counting in Chapter 3 are useful in estimating the computer resources required. A fair estimation of the time for an optimization is required before coding could be started. 5.4 Constrained Optimization of the GMR In this section, we explain how to apply simulated annealing optimizing the giant magnetoresistance. As discussed in Chapter 3, there are two definition of the giant magnetoresistance. Here we adopt the "inflated" definition of the GMR: GMR = Â— -1, (5.3) because it is "optimization friendly". If we were to use the alternate definition, GMR = 1 Â— aAp/crp, the GMR would be bounded by 100% and a significant improvement from, say, 80% and 85% would look small compared to the aimeahng temperature, and could be washed away until the temperature is very low, hence extending the annealing time significantly. When the physical property to be optimized is "optimization unfriendly," one needs to find a monotonic function of the property that magnifies the regime of interest. Not every configuration qualifies as a GMR structure. Unqualified configurations, such as those with no magnetic layers, will have zero or very small GMR. When the computational cost for each configuration is small, one can rely on the exponentially small probability to make sure the unquahfied configurations are not chosen. However, when the computational cost for each configuration is large, the time wasted in computing the unqualified structures become very costly. In com- PAGE 99 91 puting the transport properties, self-consistencies in the self-energy and the vertex correction are required, therefore the result for the last configuration is not useful in finding the result for the present configuration. Therefore each GMR calculation is expensive. One can also use a test algorithm to reject unqualified configurations, but it would still waste significant amount of time when the qualifying configuration subspace is much smaller than the full configuration space. Therefore, it is better to perform the optimization in the constrained subspace of GMR structures. The algorithm of generating the configurations in the constrained subspace will be more comphcated, but it will produce results in a much shorter time. To optimize the GMR using simulated annealing, we consider the constrained subspace of configurations consisting of a unit cell with two magnetic layers separated by two non-magnetic spacer layers. After an arbitrarily chosen initial unit cell, each subsequent configmation is generated by one of the following Monte Carlo moves: (i) inserting a monolayer, (ii) removing a monolayer, or (iii) changing the composition of a monolayer. To make sure we do not leave the constrained subspace of GMR, a move is rejected unless it will not eliminate a layer, or turn a magnetic layer into a non-magnetic layer by, for example, changing the atoms from Co to Cu. To prevent complications such as pin holes and magnetic dead layers, we further limit the minimum thickness of a layer to two mono-atomic layers. If the GMR of the new configuration is higher than the previous configuration, then it is accepted. On the other hand, if the new GMR is lower, it is only accepted with a probabihty of exp(AGMR/T), where AGMR is the change in the GMR and T is the simulated annealing temperature. The process is continued as PAGE 100 92 the annealing temperature is gradually lowered until there is no fiirther change in the configuration and a final near-global maximum in the GMR is attained. 5.5 Optimal GMR Configuration in CoNiCu superlattices In the following, we discuss the parameter and results of a study looking for the optimal configuration for superlattices made with Co, Ni, and Cu, in conditions similar to GMR experiments at the room temperature. To make the study more manageable, we only consider the fcc(lll) lattices. As mentioned in Chapter 2, the tight-binding parameters for Equation 2.1 are obtained from fits to density-function calculations where up to second-nearestneighbor hoping energies are included. In order to simulate experimental conditions, we fixed the bulk, spin-independent scattering parameter ^buik at 0.2 eV^. Using this parameter, the resistivity for bulk Cu is calculated as 3.1 /iQcm, about twice that of the room temperature resistivity of Cu. Therefore, we are modeUng a disordered metal system. The spin-dependent impurity potential caused by interface diffiision is found by fitting the tight-binding model to a supercell densityfunctional calculations of an impurity atom in the host. The values for the impiurity potentials obtained are as follows: for the majority spin, Awco/Cu =1-7 eV, AwNi/Cu =0.56 eV, Awco/Ni =0.56 eV; and for the minority spin, Attco/Cu =2.9 eV, AtiNi/cu =1-05 eV, Auco/m =0.64 eV. The interface mixing percentage is fixed in the follow way: for each monolayer, 10% of the atoms diffuse into each of the two monolayers above and below. Using the above parameters, the GMR for [C03/CU3/C03/CU3] superlattice is 36%. This value is similar to experiments. The exact value for the mixing percentage in the experiment is imknown because it is PAGE 101 93 very difficult to measure. However, as we will see, the optimal atomic configuration found in this work is not very sensitive to the errors in the parameters because the effect is robust. Figure 5.2 shows the GMR as a function of the simulated annealing step. The initial configuration (a) is [C03/CU3/C03/CU3], a symmetric structure made of four layers, each made of three monolayers of either Co or Cu. Its GMR is calculated to be 36%. However, the initial configuration does not affect the end results. Initially, the annealing temperatiue T is set to be 60%. As the optimization progresses, T decreases slowly and reducing the probability of accepting configiirations which decrease the GMR. As we can see from Figure 5.2, the GMR goes up to about 200% and then goes down to less than 50%. But at step 89, the GMR jumps to 450% and no other configiirations are accepted because any change in the configuration would have reduced the GMR too much compared with T. We continued the optimization to 600 steps and no fiurther change was observed. The atomic configuration for the highest GMR (d) is [Ni2/Cu2/Ni2/Cu2], its GMR is 450%. This is a very surprising result because most room temperature studies to date have GMR less than 100%. The GMR is sensitive to the band structure and hence the atomic configurations. For example, configurations (b) is [Ni2/Cu3/Ni2/Cu2], its GMR is 221%; configuration (c) is [Ni2/Cu2/Ni2/Coi/Cu2], its GMR is 38%. Both configurations differ from the optimal configuration (d) by only one monolayer; however, their GMR values are drastically different. PAGE 102 94 50 100 annealing step 150 Figure 5.2: GMR as a function of the simulated annealing step. An initial configuration (a), which does not affect the result, is first chosen. To optimize the GMR with simulated annealing, new configurations are generated by Monte Carlo moves while the annealing temperature T is slowly lowered. When there is no more change in the GMR, the process stops and the final configuration (d) is the optimal configuration. The four marked configurations and their GMR values are: (a) [C03/CU3/C03/CU3] 36%, (b) [Nis/Cus/Nia/Cua] 221%, (c) [Ni2/Cu2/Ni2/Coi/Cu2] 38%, and (d) [Nis/Cus/Nia/Cua] 450%. Notice that the GMR is very sensitive to the configuration. Both (b) and (c) differ from (d) by only one monolayer; however, their GMR values are very diflFerent. 5.6 Understanding the Large GMR In Chapter 4, we show that impurity atoms and structural disorder often reduce the giant magnetoresistance significantly. Therefore, it is surprising to have large GMR for similar the disorder parameters. We note that the GMR as a function in the high dimension configuration space is complicated. Therefore, it is possible that some regime of the configuration space has been overlooked. To get a better picture of the GMR as a function in the configuration space, we study the PAGE 103 95 statistical distribution of the GMR. To do this, we compute the GMR of a large number of randomly generated configurations that satisfies the criteria used in the last section. In order to obtain a reasonably accurate distribution, the number of configuration in a statistical study is a lot larger than the number of steps in the simulated annealing. We study two sets of configurations, each has 300 configiurations, their statistical distributions are similax. Therefore, we believe the number of configurations is large enough to give an reliable pictmre of the GMR distribution. On our Linux cluster, it takes about two weeks for the calculations. The histogram of the GMR of the 600 configurations is plotted in Figxure 5.3. The average GMR is 22%, and the standard deviation is 36%. The most important feature of this histogram is that there are two groups of configurations. About 97% of the configmrations are located near the main peak, they have GMR less than 100%. On the other hand, about 3% of the configurations have GMR many standard deviations away from the main peak. They scatter in the range from 100% to 450%.^ Most of the large GMR configurations are ultra-thin Ni/Cu superlattices, only two of the large GMR configurations contain Co. Since the large GMR configurations are rare, they may have been overlooked in previous studies. The very different statistical distributions in the two group suggests there is a significant physical reason behind. In order to understand the difference, in Figure 5.4 we study two typical samples in each group, [Ni2/Cu2/Ni2/Cu2] (circles) for large GMR, [C02/CU2/C02/CU2] (crosses) for regular GMR. We plotted the Figure 5.4(a) Â£?-density of states for spin-up channel of the parallel magnetic configiuration, ^The optimal configuration, which was found in the last section and has a GMR of 450%, was not among the 600 randomly generated configurations because of the large munber of possible configurations. PAGE 104 96 T 1 r ' n n r -TOO 0 100 200 300 GMR (%) Figure 5.3: Histogram of the GMR in the configuration space. In order to get a better picture of the GMR in the configuration space, we study the statistics of the GMR for 600 randomly generated configurations. The average GMR is 22%, and the standard deviation is 36%. As we can seen, about 97% of the configurations are located near the main peak. They have GMR less than 100%. On the other hand, the GMR of about 3% of the configurations are many standard deviations away from the main peak, they scatter in the range from 100% to 450%. Most of the large GMR configurations are ultra-thin Ni/Cu superlattices. Only two of the large GMR configurations contain of Co. Since the large GMR configmrations are rare, they may have been overlooked in previous studies. Figure 5.4(b) the conductivity for the same channel, and Figiu-e 5.4(c) the GMR, as fmictions of the Fermi energy in the rigid band picture. The physical Fermi energy corresponds to Â£'f = 0 Ry. As seen in Figure 5.4(a), the two configurations have similar band structures. The main difference is that while [Ni2/Cu2/Ni2/Cu2] has a sharp d-DOS despite of the 10% interface mixing and 0.2 eV^ bulk scattering, [C02/CU2/C02/CU2] has a spread d-DOS even at 5% interface mixing and the 0.2 eV^ bulk scattering. This is because the scattering between spin-up Co and Cu ^ 80 q *--' CO .gj 60 c o o 0 40 a> .0 E 1 20 PAGE 105 97 is about a three times the scattering between spin-up Ni and Cu. Therefore, the scattering in Co/Cu interface per atom is about 9 times as effective as the scattering in the Ni/Cu. As a result, [C02/CU2/C02/CU2] has a spread d-peak and its (i-DOS at the Fermi energy is large, while [Ni2/Cu2/Ni2/Cu2] has a sharp c?-peak and its d-DOS at the Fermi energy is small. As seen from figiures 5.4(a) and 5.4(b), the conductivity increases as the c?-DOS decreases. This is because when d-DOS is large, s-d hybridization reduces the average Fermi velocity significantly, and hence reducing the conductivity. Since the conductivities in the spin-down channel and in the anti-parallel magnetic configiu-ation do not vary as much as this channel, the GMR curve roughly follows the conductivity curve for the spin-up channel of the parallel magnetic configuration. We notice that both examples eventually have large GMR at high enough energy, however, only the GMR at the Fermi energy, u; = 0, is measiured in many experiments. As a result, there is a qualitative difference in the GMR of [Ni2/Cu2/Ni2/Cu2] and [C02/CU2/C02/CU2] despite of their similar band structure. ^ 5.7 Conclusion Using a highly optimized GMR code, we have used the simulated annealing algorithm to search for the optimal atomic configuration for CPP-GMR. For conditions similar to experiments at the room temperatmre, we foimd that [Ni2/Cu2/Ni2/Cu2] has the highest GMR, which is 450%. This is smrprising because GMR measurements for most configurations are usually much lower when scattering is large. We have also studied the statistical distribution of GMR, which has given us insight about GMR as a fimction in the configuration space. The configurations are found PAGE 106 98 >; 40 EE o b 5 0.05 t*e 0 1000 5 500 a -8: n r ^ Â— I Â— h ! M ' 1 * Â«Â» * 0 ? 04 -0.03 -0.02 -0.01 0 0.01 0.02 (o(Ry) Figure 5.4: Comparison of [Ni2/Cu2/Ni2/Cu2] (circles) and [C02/CU2/C02/CU2] (crosses) superlattices. Although their band structure is similar, the NiCu superlattice has a large GMR (450%), while the CoCu superlattice has a smaller GMR (50%). This difference is caused by the spreading of the rf-peak from lower energies to the Fermi level, u; = 0. (a) The d-density of states of the spin-up channel for the parallel magnetic configuration, c?-DOS||, of the [Ni2/Cu2/Ni2/Cu2] at the Fermi level is much lower then that of the [C02/CU2/C02/CU2]. (b) This leads to a higher conductivity, a||, for the same spin channel of [Ni2/Cu2/Ni2/Cu2] than in [C02/CU2/C02/CU2]. Since the conductivities a^, a-^, cr|| are relatively fiat compared with an near the Fermi energy, the GMR (c) resembles cr-fi (b). to be in two groups, 97% of the configurations have regular GMR, and 3% of the configurations have large GMR. Since the large GMR configurations are rare, they might have been overlooked. We also show that the large GMR is caused by the relatively sharp d-peak and hence low o?-DOS at the Fermi energy for the spin-up channel of the parallel configuration, which in turn causes high conductivity in the same channel and a large GMR. PAGE 107 99 We have studied more than a thousand configurations, however, the cost of human resource is low. Once we learned that large GMR exists in a certain regime, it took us Uttle time to study this regime and find out the physics behind. In a sense, the computer decides what to look for, in a tedious but slightly 'intelligent' way, yet it is effective in terms of actual cost for the research project. It would be a waste, if at all possible, to have people performing a search the way the computer did. To conclude, the optimization algorithms has made our research more effective. The effect found here is a robust effect with simple physical explanation, yet we would not have found it with conventional discovery process. PAGE 108 CHAPTER 6 CONCLUSION In this final chapter, we present a brief review of the dissertation and discuss its implications as well as future directions. In Chapter 2, we explained how to obtain the parameters for the multi-band tight-binding Hamiltonian and for the disorder. We also discussed how to calculate the conductivity within the model, including the vertex correction necessary for current conservation. In Chapter 3, we developed a highly optimized parallel code capable of finding the self-energy in 0{N) and the conductivity in 0{N'^). In Chapter 4, we discussed the some of the experiments in the current-perpendicular-to-plane giant magnetoresistance and applied oiu: method to configurations similar to the experiments. Our calculations agrees with the resistivity and giant magnetoresistance data of Cyrille et al, and we can explain quantitatively data of Gijs et al. by a shunt resistor model. Although we only studied the magnetic multilayer superlattices, our method can be applied in a wide variety of systems and more complicated geometries. With small modification, this code should also work for non-periodic systems, such as systems with conducting leads or systems with complicated interface structures. Most importantly, the size dependence will remain the same. Also, the fast conductivity code may also be used to study large systems. In Chapter 5, we turned to a new direction in our research. Using the simulated annealing method, we searched for the optimal giant magnetoresistance structiues made with Co, Ni, and Cu monolayers. We found that some ultra-thin Ni/Cu su100 PAGE 109 101 perlattices have a high GMR, as much as 450%, even when the bulk and interface disorder are moderate. This is a surprising discovery because it is expected that disorder at this level should reduce the GMR significantly. The large GMR can be understand by noting that in a small number of ultra-thin Ni/Cu superlattice, it is possible to have a low majority d-density of states at the Fermi level even with moderate disorder. In 3d transition metals, low d-density of states at the Fermi level implies low resistivity in this channel, and hence a large giant magnetoresistance. The effect found in this study is not very sensitive to the parameters because it is robust. Although this effect is large and has a clear physical explanation behind, it was not discovered before because the regime of this effect is very small compared with the size of the configuration space. With the simulated annealing, the computer is "intelligent" enough to avoid trapping by local optima while using the local information to find the global optima. This method makes our research more efficient and eventually leads to findings that are not likely to be discovered with conventional methods. As far as we know, this is the first optimization study of transport property at the atomic level. By shifting some tedious task to the fast performing computers of today, one can make a qualitative difference in the process of scientific discovery. This is an important trend, although what we used is very far away from the "engines of design" predicted in Eric Drexler's Engines of Creation [66], which "will be able to create bold new design without human help". We still need to ask the right question, formulate the physical model, verify, understand, and make use of the PAGE 110 102 results. Nevertheless, a significant amount of hiunan resources which would have spent in trial-and-error is reduced. The program used here may also be apphed to optimize the GMR for hot electrons, where the measurement is made at energies higher than the Fermi energy. This approach may also be applied to the area of spin transport across semiconducting/metaUic interfaces, where the band mismatch is very large and optimal spin transport structmes may not be find intuitively. With a suitable model Hamiltonian for the semiconductors, it may be possible to optimize spin transports across semiconducting/metaUic interfaces. We expect more types of problems can be studied with this approach as faster computers are available. PAGE 111 APPENDIX A FAST INVERSE ALGORITHMS To compute the Green's function efficiently, it is essential to have fast inverse algorithms that converge well. If the size of a tight-binding Hamiltonian is A'^, we would like to make use of the sparsity of the obtain the full Green's function in 0{N'^). In addition, since only the onsite Green's function are needed to find the self-energy, it is desirable to have an algorithm that gives the diagonal of the inverse in 0{N). The matrices of interest are semi-hermitian^ tridiagonal block matrices, with or without periodic boundary conditions. Although the inverse of such matrices is common in physics applications, there are not many inverse algorithms or codes targeting this type of matrix, especially with periodic boundary conditions. In Section A.l, we review the algorithm discribe by E. M. Godfrin in 1991 to find the inverse of block tridiagonal matrices, without periodic boundary conditions. It costs^ 0{N'^) to compute the full inverse and 0{N) to compute the diagonal of the inverse. ^Matrix M is a semi-hermitian matrix if M{j,j') = M{j',j)* for all j ^ j'. For example, M = u; Â— H + iri is a, semi-hermitian matrix. ^Since the elements in block matrices are themselves matrices, element additions are much cheaper than multiplications, and are hence neglected in cost counting. This method of counting operations is different from that of general matrices, where both a multiplication and an addition are count as floating point operations(flops) and cost the same in many computers with mathcoprocessors. In this chapter, the elements are m x m matrices, and the operation cost counts are therefore in unit of flops. 103 PAGE 112 104 In Section A. 2, we develop a method to solve the inverse of tridiagonal block matrices with periodic boundary conditions. In many inverse algorithms such as the Godfrin algorithm, it is crucial to have fixed boundaries from which iterations start. Although the periodic boundary condition only adds two blocks to the original matrix, it makes finding the inverse much more difficult. To overcome this difficulty, we use the Woodbury formula to compute the change in the inverse of a matrix after a rank one change. Using both Godfrin and Woodbury methods, we are able to compute the inverse of a block tridiagonal matrix with periodic boundary condition and the diagonal of the inverse in 0(N'^) and 0{N), respectively. In Section A. 3, we improved Godfrin's method for inverse for both types of matrices. Detailed operation counts are provided to show the importance of the results presented in this chapter. The methods in this appendix are tested against Matlab for large matrix inverses to double precision accuracy. For the type of matrices discussed here, these algorithms should converge better than the algorithms for general matrices. A.l Godfrin Inverse Algorithm for Tridiagonal Block Matrix Godfrin's method [67] is similar the one derived by Economou-Cohen, but Godfrin's result is exact and converges very well, and the formulation is particularly suitable for our calculations. The original formulation is for a semi-hermitian matrix, however, it can be formulated for a general tridiagonal block matrix. Let M be a block tridiagonal matrix, and G = be the inverse. It is best to demostrate how to compute the inverse with a 4 x 4 matrix, and then present PAGE 113 105 the general algorithm. The explicitly definition of the inverse is given by \ / Gi4 \ / 0 0 \ Mu Mu 0 0 Gn Gn Giz 1 0 M21 M22 M23 0 G21 G22 G23 G24 0 1 0 0 0 M32 M33 M34 Gsi Gz2 G33 G34 0 0 1 0 V 0 0 M43 M44 / V G41 G42 G43 G44 / V 0 0 0 1 / This is a set of 16 equations; however, the equations for each column of G can be solved separately. For example, the equations for column 2 are: M11G12 + M12G22 = 0 (A.l) M21G12 + M22G22 + M23G32 = 1 (A.2) M32G22 + M33G32 + M34G42 = 0 (A.3) M43G32 + M44G42 = 0. (A.4) This set of equations can be solved independent of other equations. From Equation A.4, G42 = -M4^^M43G32 M34G42 = -M34M4-/M43G32 = -X3G32. (A.5) PAGE 114 106 Substituting the result into Equation A. 3, we have [Ms3 Â— XzlGsi = Â— M32G22 ^23^32 = -M23[M33 -X3]-^M32G'22 = -X2G22. (A.6) We leave behind Equation A. 2, which is non-zero on the right hand side, until other equations are used. The idea is that we first represent every variables above and below G22 in terms G22, then Equation A. 2 can be solved easily. Prom Equation A.l, we have G12 = -Mfi^Mi2G22 . . . M21G12 = -M2iMf/Mi2Gi2 = -I2G22. ^ : / , ' (A.7) Substitute the Equations A.6 and A. 7 into Equation A. 2, we obtain Â— ^2^22 + M22G22 Â— X2G22 = 1 ' G22 = [M22 X2 Fa]"'(A.8) The rest of the variables are: G42 = -M^M4iGz2 PAGE 115 107 Gn = -Mf/Mi2Gi2. (A.9) We have spent 0{N) operations in order to solve this set of A'^ equations. However, since the form of the equations on the left side of Equations A.1-A.4 are the same for other columns, a large part of the result can be reused when solving equations for column 1, 3 and 4. The formal procediu-e is as follows : Define : Xn = 0 Xi = Mi,i+i[Mj+i,i+i Xi+i]-^Mi+i^i for 1 < i < 1 n = 0 Yi = Mi,i_i[Mi_i,i_i Fi_i]-iMi_i,i for 2 < i < Ci = -[Ma Xi]-'Mi^i_i for 2 < i < A = -[Mii Fi]-'Mi,i+i for 1 < i < A^ 1, (A.IO) then Gii = [Mu-Xi-Yi]-' Gij = CiGi-ij for i > j Gij = Did+ij. for i PAGE 116 108 The diagonal and any column j of Gij can be found in 0{N).^ The cost of the full inverse is A^^ + 0{N). A. 2 Inverse of Tridiagonal Block Matrix with Periodic Boundary Condition Computing the inverse of tridiagonal block matrix with periodic boundary condition is more complicated. We are still able to find the inverse in 0{N^) and the diagonal of the inverse in 0{N). However, the coefficients will be larger. In the following, we first discuss the Woodbury formula, which is used to compute the change to the inverse of a matrix after a rank one change. A. 2.1 Rank-one changes: Woodbury formula Let M, B he N X N block matrices, u, v be size N column block vector, the Woodbury formula [68] ^ states that if M = B-uv'^, then M-^ = B-^ + B-\av'^B-\ (A.12) where a = (1 Â— v'^B~^u)~^ has the same dimension of an element block. To prove this formula, one considers the product {B-uv'^){B-^-^B-'uaii^B-^) ^Note that in order to find a row in 0{N), we must modify the algorithm, as discussed in A.2.2. ^It is called the ShermanMorrison formula when the size of a block is 1 x 1. PAGE 117 109 uv ''^B-^uav'^B= 1 Â— uv <'^B-^ +u{l-v <'^B-^u)av'^B -1 % ... (A.13) because {l-v'^B~^u)a = 1. Similarily, the product {B~^+B~'^uav'^B''^){B-uv'^) is also equal to 1. Although this formula looks like a perturbation method, it is indeed exact. A.2.2 Applying the Woodbury Formula The Godfrin algorithm in A.l rehes on fixed boundary conditions in order to start the iteration. For a.NxN block tridiagonal matrix M with periodic boundary conditions, there are two more block elements: Mijv and Mat i . This two elements makes the inverse much harder to solve than it looks because there is no place to starting the iteration. In order to find the inverse, we first write M into a tridiagonal block matrix without the boundary condition plus the periodic conditions. The inverse of the first term is computed by the Godfrin method, and the correction due to the periodic condition is computed by the Woodbm-y formula. In other words, write M in the following form: . . M = B-uv'^ PAGE 118 110 Mu MiN Mi2 ... 0 M21 M22 0 0 0 \ 0 0 -MiN 0 0 -Mm 0 0 1 0 MiV-l,Ar-l Mn~i,n Mn,n-i Mnn Mni j 0 1 (A.14) where B is the tridiagonal matrix without periodic boundary condition, u is the column vector, and is the row vector. B'^ can be found by Godfrin's algorithm discussed in section A.l. can then be computed by the Woodbury formula discussed in A.2.1. A. 2. 3 Cost Considerations The method in A. 2. 2 looks very simple, however, our goal of computing the diagonal of the inverse in 0{N) can be acheived only by examining the method carefully. Since u and v in Equation A.14 are sparse, we should write down the Woodbury formula explicitly for our case. In the equation a = {l-vlBr^\)-\ (A.15) PAGE 119 Ill where repeated indices are summed automatically, we only need to calculate {B^,^ + B]:j\)ui + {B^j^UN + B]^l!)uN, (A.16) It costs only 6 multiplications. The correction term in Equation A. 10 is more complicated: {B-'uav''B-% = {B-'uaUv^B-% (A.17) where {B~'^ua)i = Bii{uia) + BiNiuNOt), {v^B-% = vJ{B-%+vl{B-')^j = {B-% + {B-')r,j. (A.18) If we do this carefully, i.e. multiply un and ui by q first, it costs AN + p multiplication, where p is the number of elements we want. Therefore, the cost for the correction, {B~^uav^B~^)ii, to the diagonal is 5N, and the cost for the correction to the whole matrix is N'^ + AN. In Equation A.16, we need both columns {B''^)ii and (J5~^)ijv, and rows {B~^)ij and {B~^)f^ij. However, in section A.l, we only know how to obtain {B~^)ij, where j varies, in a row form. If we want to obtain {B~^)ii^, then for each i, it takes 0{N) steps, and the whole column costs 0{N'^). In order to obtain the columns PAGE 120 112 also in 0{N), we need to find the inverse of the tridiagnonal matrix in the cohimn form. Similar to Equations A. 10 and A.ll, we have q = -Mi_,,i[Mu-X, 1-1 for 2 j. Note that Cj/ (Cj)^ because Ma is usually semi-hermitian or worse. (A.20) A. 3 Efficient Coding for the Godfrin Algorithm We can speed up the original formula listed in Godfrin's paper by being careful not to repeat some of the calculations. Also, both the column and the row Godfrin inverse can be computed at the same time. The proper procedure is as follows: Xn = 0, Wr, = [Xn-Mmn]-\ Cn = WnMn,n-u C'j^ = Mn-i,nWn PAGE 121 113 and Y, = 0, D, = Z,Mn, D[ = M21Z1 Y2 = M2i[Mn-Yi]-'Mu = -M21D, = -D[Mn. (A.21) To calculate Xj,Yj,Wj,Zj,Cj, Dj for all j, it takes 6A'' Â— 2 multiplication or inverse operations. It takes A'^ more inverses to find out the diagonals. Therefore, the cost for finding the inverse of a tridiagonal block matrix by Godfrin algorithm is p + 67V + 0(1), where p is the number of elements we need. It takes 2N more calculations to find all the C'^Dj, and 4A'^ Â— 2 operations to find out Gij,GNj,Gii,GiNCounting the cost of both the Godfrin and the Woodbury methods, the total cost to calculate p elements of the inverse of a tridiagonal block matrix with periodic boundary condition is 2p + 16N + 0(1). PAGE 122 114 This method is very efficient in finding the diagonals when TV is larger than about 18. Since these methods are designed specifically to solve the tridiagonal block matrices, it should be more efficient and accurate than Gaussian elimination with or without pivoting. To find the inverse for tridiagonal block matrix with periodic boimdary condition and 9x9 element blocks, it takes the Matlab sparse matrix inverse rountine roughly 9.4A^^ x 9^ flops, as opposed to roughly 2N'^ x 9"^ flops in this method. PAGE 123 APPENDIX B COHERENT POTENTIAL APPROXIMATION The coherent potenetial approximation (CPA) [26, 49-50] is an average method for computing the density of states in substitution alloys within the singleimpurity approximation. Two commonly quoted approaches are the effective medium approach and the locator approach, which lead to the same set of equations. Consider an atom of type j in an effective medium with potential VeIn real space matrix notation, the Green's function for the effective medium is where we have written the Hamiltonian as a sum of the offsite part i/offsitej and the onsite part 14 which is adjusted later. The T-matrix due to the atom is Go=[E //offsite -Ve + IT]] -1 (B.l) (B.2) (B.3) where aj ^ [7 {V, K)Go] -1 (B.4) 115 PAGE 124 116 If we consider only a single-site impurity ^ and onsite impurity potentials, then Tj is non-zero only at a single site. This simplifies the problem significantly. Now, instead of solving a problem for the entire disordered bulk system, we only have to deal with a single-site T-matrix (a single matrix block/element in the full matrix). The effect on the whole medium will be carried out by adjusting the effective potential and re-applying Equation B.l. w''-^ The CPA requires the ensemble average scattering to be zero at each energy E, i.e. This condition can be satisfied by modifying Vg. The new Ve can be found by rewriting the above equation as For a multiband Hamiltonian, the order of operations can not be changed because these are matrix equations. Equations B.l and B.6 are then applied repeatedly to obtain a selfconsistent solution for Ve and Gq. The above algorithm is tested to work for a 9-band tightbinding model of fee NiFe and bcc FeNi. The effective potential K looks like a weighted average of Vi and V2 where the energy dependent weight is adjusted by Go and K itself. Obviously, the above equation works well in the pure host limits, Ci = 0 or C2 = 0. ^Thus CPA neglects cluster effects. (T) = ciTi + C2T2 = c,ai{V, K) + C2Â«2(^2 K) = 0. (B.5) Ve = [Cifti -I020:2] ^{ciCliVi + 0202^2). (B.6) PAGE 125 117 After obtaining the effective potential and Green's function of the effective medium, we can compute the Green's function for a specific sites. For an atom of type j, the Green's fimction is Gj = Go + GoTjGo. (B.7) This equation is important for finding the site-specific quantities such as the local density. PAGE 126 REFERENCES [1] G. Prinz. Physics Today, 48(4) :58, 1995. [2] F. J. Himpsel, J. E. Ortega, G. J. Mankey, and R. F. Willis. Adv. Phys., 47(4):511, 1998. [3] M. N. Baibich, J. M. Broto, A. Fert, F. Nguyen Van Dau, F. Petroff, P. Etienne, G. Creuzet, A. Friederich, and J. Chazelas. Phys. Rev. Lett, 61:2472, 1988. [4] M. A. M. Gijs, S. K. J. Lenczowski, and J. B. Giesbers. Phys. Rev. Lett, 70:3343, 1993. [5] M. A. M. Gijs, J. B. Giesbers, M. T. Johnson, J. B. F. aan de Stegge, H. H. J. M. Janssen, S. K. J. Lenczowski, R. J. M. van de Veerdonk, and W. J. M. de Jonge. J. Appl. Phys., 75:6709, 1994. [6] M. C. Cyrille, S. Kim, M. E. Gomez, J. Santamaria, K. M. Krishnan, and I. K. SchuUer. Phys. Rev. B, 62:3361, 2000. [7] M. C. Cyrille, S. Kim, M. E. Gomez, J. Santamaria, C. Leighton, K. M. Krishnan, and I. K. SchuUer. Phys. Rev. B, 62:15079, 2000. [8] W. Varra. Appl. Phys. Lett, 66:2579, 1995. [9] W. P. Pratt Jr, S. F. Lee, J. M. Slaughter, R. Loloee, P. A. Schroeder, and J. Bass. Phys. Rev. Lett, 66:3060, 1991. [10] W. P. Pratt Jr, S. F. Lee, Q. Yang, D. Holody, R. Loloee, R A. Schroeder, and J. Bass. J. Appl. Phys., 73:5326, 1993. [11] Q. Yang, S.-F. Lee, P. Holody, R. Loloee, P. A. Schroeder, W. P. Pratt Jr, and J. Bass. Physica B, 194-196:327, 1994. [12] Q. Yang, P. Holody, R. Loloee, L. L. Hemy, W. P. Pratt Jr, P. A. Schroeder, J. Bass, and P. A. Schroeder. J. Mag. Magn. Matter., 126:406, 1993. [13] T. Ono, Y. Sugita, K. S. K. Mibu, N. Hosoito, and T. Shinjo. Phys. Rev. B, 55:14457, 1997. 118 PAGE 127 119 [14] W. Oepts, M. A. M. Gijs, A. Reinders, R. M. Jmgblut, R. M. J. van Gansewinkel, and W. J. M. de Jonge. Phys. Rev. B, 53:14024, 1996. [15] L. Piraux, J. M. George, J. F. Despres, C. Leroy, E. Ferain, R. Legras, K. Ounadjela, and A. Fert. Appl. Phys. Lett, 65:2486, 1994. [16] T. Valet and A. Fert. Phys. Rev. B, 48:7099, 1993. [17] S. Zhang and P. M. Levy. J. Appl. Phys., 69:4786, 1991. [18] Z.-G. Zhang and W. H. Butler. J. Appl. Phys., 81:4576, 1997. [19] S. Zhang and P. M. Levy. Phys. Rev. B, 57:5336, 1998. [20] J. Barnas and A. Fert. Phys. Rev. B, 49:12835, 1994. [21] A. Vedyayev, N. Ryzhanova, B. Dieny, P. Dauguet, P. Gandit, and J. Chaussy. Phys. Rev. B, 55:3728, 1997. [22] E. Yu. Tysmbal and D. G. Pettifor. Phys. Rev. B, 54:15314, 1996. [23] J. Matton, A. Umerski, and Villeret. Phys. Rev B, 55:14378, 1997. [24] S. Krompiewski, M. Zwierzycki, and U. Krey. J. Phys.: Condens. Matter, 9:7135, 1997. [25] J. Chen, T.-S. Choy, and S. Hershfield. J. Appl. Phys., 85:4551, 1999. [26] T.-S. Choy, J. Chen, and S. Hershfield. J. Appl. Phys., 86:562, 1999. [27] K. M. Schep, J. B. A. N. van Hoof, P. J. Kelly, G. E. Bauer, and J. E. Inglesfield. Phys. Rev. B, 56:10805, 1997. [28] Q. Yang, P. Holody, S.-F. Lee, L. L. Henry, R. Loloee, P. A. Schroeder, W. P. Pratt Jr, and J. Bass. Phys. Rev. Lett., 72:3274, 1994. [29] J. Bass, Q. Yang, S.-F. Lee, P. Holody, R. Loloee, P. A. Schroeder, and W. P. Pratt Jr. J. Appl. Phys., 75:6699, 1994. [30] T.-S. Choy, J. Naset, J. Chen, S. Hershfield, and C. Stanton. Fermi surface database, http://www.phys.ufl.edu/fermisurface, accessed 7/1/2001. [31] T.-S. Choy, J. Naset, J. Chen, S. Hershfield, and C. Stanton. Bull. Amer. Phys. Soc., Mar. Meet, 45(1):L36 42, 2000. PAGE 128 120 [32] F. Darema, S. Kirkpatrick, and V. A. Norton. IBM Journal of Research and Development, 31:391, 1987. [33] M. Mitchell. An Introduction to Genetic Algorithms. MIT, Boston, 1995. [34] J. Koza. Genetic Programming. MIT, Boston, 1992. [35] J. R. Chelikowsky. J. Phys. D: Appl. Phys., 33:R33, 2000. [36] D. M. Deaven and K. M. Ho. Phys. Rev. Lett, 75:288, 1995. [37] A. Franceschetti and A. Zunger. Nature, 402:60, 1999. [38] R. Kubo. Can. J. Phys., 34:1274, 1956. [39] R. Landauer. IBM J. Res. Dev., 1:223, 1957. [40] R. Landauer. Ma^., 21:863, 1970. [41] W. Kohn and L. Sham. Phys. Rev. A, 140:1133, 1965. [42] S. Ohnishi, A. J. Freeman, and M. Weinert. Phys. Rev. B, 28:6741, 1983. [43] G. M. Stock, Yang Wang, , D. M. C. Nicholson, W. A. Shelton, W. M. Temmerman, Z. Szotek, B. N. Harmon, and V. P. Antropov. Materials Theory, Simulations, and Parallel Algorithms. Symposium, xiv+611:157, 1996. [44] K. Majumdar. Study of Transport Properties in Magnetic Nano structures. Ph.D. Dissertation, University of Florida, Gainesville, 1999. [45] Rickayzen. Green's Functions for Condensed Matter Physics. Academic, NewYork, 1984. [46] S. Hershfield. Ann. Phys., 196:1, 1989. [47] G. D. Mahan. Many-Particle Physics. Kluwer Acadenaic/Plenum Publishers, New York, 2000. [48] D. A. Papaconstantopoulos. Handbook of the Band Structure of Elemental Solids. Pleunum, New York and London, 1986. [49] R. J. Elliott, J. A. Krumhansl, and P. L. Leath. Rev. Mod. Phys., 46:3, 1974. [50] F. Yonezawa and K. Morigaki. Suppl. Prog. Theor. Phys., 53:1, 1973. PAGE 129 121 J. J. Dongarra, J. Du Croz, I. S. Duff, and S. Hammarling. ACM Tran. Math. Soft., 16:1, 1990. E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK User's Guide. SIAM, Philadelphia, 1992. Bjarne Stroustrup. The C++ Programming Language. AddisonWesley, Reading, 1995. R. Pozo. Journal of Supercomputing Applications and High Performance Computing, 11.3, 1997. R. Pozo. Template numerical toolkits, http://math.nist.gov/tnt/index.html, accessed 7/1/2001. T. L. Veldhuizen. Proceedings of the 2nd International Scientific Computing in Object Oriented Parallel Environments, xi-|-242:223, 1998. T. L. Veldhuizen. Blitz-F-f. http://www.oonumerics.org/blitz/, accessed 8/5/2001. S. S. Parkin, R. Bhadra, and K. P. Roche. Phys. Rev. Lett, 71:1641, 1991. S. S. Parkin. Phys. Rev. Lett, 67:3598, 1991. P. Griinberg, R. Schreiber, Y. Pang, M. B. Brodsky, and H. Sowers. Phys. Rev. Lett, 57:2442, 1986. W. Schepper, K. Diplas, and G. Reiss. J. Appl. Phys., 87:6597, 2000. M. Freyss, D. Stoeffler, and H. Dreysse. Phys. Rev. B, 54:18, 1996. C. Kittel and in. Solid State Physics, 22:1, 1968. L. M. Roth, H. J. Zeiger, and T. A. Kaplan. Phys. Rev., 149:519, 1966. S. R. Wilson and A. W. Czarnik. Combintorial Chemistry: Synthesis and Application. Wiley, New York, 1997. K. E. Drexler. Engines of Creation. Anchor, New York, 1986. E. M. Godfrin. J. Phys. : Condens. Matter, 3:7843, 1991. W. H. Hager. Applied Numerical Linear Algebra. Pentice Hall, Englewood Cliffs, 1988. PAGE 130 BIOGRAPHICAL SKETCH Tat-Sang Choy was born on December 3, 1971, in GuangZhou, China. He moved to Hong Kong in 1981. He attended the Sing Yin Secondary School, where he spent most of the time in the astronomy club and writing computer programs. He then entered the Chinese University of Hong Kong and graduated with a B.Sc. degree in physics in 1995. He was awarded the the Gold Prize of the Undergraduate Projects Awards 1995, for his work in nonlinear response in clusters with Professor Kin-Wah Yu. In the same year, he moved to the United States to study in the Physics Department of the University of Florida. In the summer of 1996, he joined Professor Selman Hershfield's research group and worked on the electron transport theory of magnetic nanostructures. He will start his postdoc at the Ohio State University this September. 122 PAGE 131 123 General Audience Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy ELECTRON TRANSPORT THEORY IN MAGNETIC NANOSTRUCTURES By Tat-Sang Choy August 2001 Chairman: Selman P. Hershfield Major Department: Physics Magnetic nanostructures are magnetic devices which have at least one dimension on the order of one nanometer, 10"^ m, and they have become the focus of a new trend in nanoelectronics because of their applications in magnetic sensors, magnetic memories, and magnetic reading heads in hard disks drives. For example, a reading head in a magnetic hard-disk is a nanostructure made by growing a magnetic layer on a substrate, followed by a non-magnetic layer, and then another magnetic layer. The thickness of these layers are of the order of one nanometer. The resistance of this structure changes as a magnetic field is applied, therefore it is used as in a reading head. In this dissertation, we first develop a method to model the electrical response of the magnetic nanostructure in the atomic scale. This method is then improved such that it can be applied in very complicated structures. Then we apply this method to study the giant magnetoresistance effect and successfully explain some of the existing experiments. Finally, we turn to a new direction of research in nanostructures in which we use an optimization method on high performance computers PAGE 132 to search for new nanostructures. We found a magnetic nanostructure that has a surprisingly large giant magnetoresistance ratio. It is hoped that this nanostructure will be realized in experiments. Using computers to search for thousands of nanostructures automatically, the efficiency of scientific discovery is increased and a significant amoimt of human resources is saved. PAGE 133 I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quaUty, as a dissertation for the degree of Doctor oj Philosophy. I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosoph;^. Arthur F. Hebard Professor of Physics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Mark W. Meisel Professor of Physics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quahty, as a dissertation for the degree of Doctor of Philosophy. Peter J. Hirsihfeld A Professor of Physics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. PAGE 134 Stephen Professor of Materials Science and Engineering This dissertation was submitted to the Graduate Faculty of the Department of Physics in the College of Liberal Arts and Sciences and to the Graduate School and was accepted as partial fulfillment of the requirments for the degree of Doctor of Philosophy. August 2001 f Dean, Graduate School xml version 1.0 encoding UTF-8 REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID EJKSKQJ0O_CNFXL2 INGEST_TIME 2015-04-23T21:27:04Z PACKAGE AA00029942_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES |