Quantum Mechanical Simulation and X-Ray Scattering Applied to Pressure-Induced Invar Anomaly in Magnetic Iron Alloy Citation Winterrose, Michael Lon (2011) Quantum Mechanical Simulation and X-Ray Scattering Applied to Pressure-Induced Invar Anomaly in Magnetic Iron Alloy. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/2CBN-7H49. https://resolver.caltech.edu/CaltechTHESIS:08302010-215155302 Abstract The Invar effect has remained at the forefront of materials research since Charles-Édouard Guillaume discovered the vanishing thermal expansion of Fe-Ni alloys in 1897. More recently, a pressure-induced Invar effect was discovered in Fe-Ni alloys, and the relationship between classical and pressure induced Invar phenomena has added complexity to the century-old struggle to comprehend the microscopic origins of Invar behavior. In this thesis I present our recent discovery of pressure-induced Invar behavior in Pd₃Fe with the ordered L1₂ structure. Nuclear forward scattering measurements show that the ferromagnetic ground state in Pd₃Fe is destabilized with pressure, collapsing around 10GPa (V/V0=0.96) to a lowspin magnetic state. From high-pressure synchrotron x-ray diffraction measurements we find a large volume collapse at ambient temperature to accompany the collapse of ferromagnetism. After the volume collapse there is a significant increase in the bulk modulus. Using nuclear resonant inelastic x-ray scattering to study the 57 Fe phonon partial density of states (PDOS) at high pressures, we find the pressure-induced magnetic transition to cause an anomalous relative softening of the average phonon frequency. Heating our sample to 650K in a furnace at a pressure of 7GPa, synchrotron xray diffraction measurements reveal negligible thermal expansion from 300 to 523 K, demonstrating pressure-induced Invar behavior in Pd₃Fe. Density functional theory calculations identify a ferromagnetic ground state in Pd₃Fe with large moments at the Fe sites. These calculations show that the application of pressure counteracts the band-filling effect of Pd. By tuning the position of the top of the 3d band with respect to the Fermi level, pressure-induced Invar behavior resembles classical Invar behavior that is controlled by chemical composition. This insight marks the first step towards a unification of our understanding vii of classical and pressure-induced Invar behavior. Pressure drives the majority-spin t 2g antibonding electronic states closer to the Fermi level. The transition to the low-spin state occurs as these t2g states move across the Fermi level, transferring charge to the minority-spin eg nonbonding electronic states. This charge transfer reduces the internal electronic pressure in the material, giving a volume reduction in the low-spin state. The movement of the t 2g states with increasing pressure results in a greater number of states at the Fermi level, increasing screening efficiency and softening the first nearest-neighbor Fe-Pd longitudinal force constants in the low-spin state. The measured and calculated magnetic transition pressures differ significantly, despite sharing similar elastic properties in both the ferromagnetic and low-spin states. The magnitude of the disagreement between theoretical and experimental magnetic transition pressures suggests a spin-disordered state exists at high pressures in Pd₃Fe. A shape discrepancy between the calculated and measured high-pressure Fe PDOS suggests significant short-range spin correlations exist in this spin-disordered state. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Invar, high pressure, diamond anvil cell, density functional theory, magnetism Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Materials Science Minor Option: Applied Physics Thesis Availability: Public (worldwide access) Research Advisor(s): Fultz, Brent T. Thesis Committee: Fultz, Brent T. (chair) Johnson, William Lewis Ravichandran, Guruswami Daraio, Chiara Jackson, Jennifer M. Defense Date: 23 August 2010 Record Number: CaltechTHESIS:08302010-215155302 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:08302010-215155302 DOI: 10.7907/2CBN-7H49 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 6017 Collection: CaltechTHESIS Deposited By: Michael Winterrose Deposited On: 16 Sep 2010 21:02 Last Modified: 09 Oct 2019 17:05 Thesis Files Preview PDF - Final Version See Usage Policy. 8MB Repository Staff Only: item control page

Quantum Mechanical Simulation and X-Ray Scattering Applied to Pressure-Induced Invar Anomaly in Magnetic Iron Alloy Thesis by Michael L. Winterrose In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 2011 (Defended August 23, 2010) ii c 2011 Michael L. Winterrose All Rights Reserved iii To Jaime, the love of my life. To my parents, Mike and Treo, my unwavering support. iv Acknowledgements This thesis was made possible through the support, collaboration, and advice of my advisor, Brent Fultz, who believed in me at a time when I did not believe in myself and set me on a path that has led me to this day. The high-pressure experimental work in this thesis was completed largely through the near superhuman efforts of Itzhak Halevy. His extraordinary effort, guidance, patience, and talents will not be forgotten. I would like to thank Channing Ahn for getting me away from my computer during the periods between the synchrotron trips and imparting his invaluable insights into the world of experimental materials work. His support and advice is much appreciated. The expert machining work of Mike Vondrus made possible the discovery around which this thesis revolves. Also, the dedicated work of the administrative assistants made our work and travel possible. I would particularly like to thank Pam Albertson and Taura Scott for their help. Many kind thanks to my committee, Brent Fultz, William Johnson, Chiara Daraio, Guruswami Ravichandran, and Jennifer Jackson. I would also like to acknowledge my undergraduate advisor, Professor J. Thomas Dickinson, for setting me on the road of computational materials physics that now finds its culmination in this thesis. For their deep knowledge and willingness to share it, their impeccable manners, and for the many interesting conversations we have engaged in over the years, I thank my Caltech office mates Rebecca Stevens, Matt Lucas, Jorge Munoz, and Chen Li. To former and current Fultz group members, Olivier Delaire, Alan Yue, Ryan Monson, Hongjin v Tan, Nick Stadie, David Abrecht, Hillary Smith, Kunwoo Kim (for showing me the funny side of hernias), Sally June Tracy, Brandon Keith, Alex Dementsov, Nikolay Markovskiy, Xiaoli Tang, Ryan Monson, Mike McKerns, Jiao Lin, Alexander Papandrew, Tabitha Swan-Wood, Shu Maio, Max Kresch, and Justin Purewal, I thank you for making my time at Caltech fruitful and enjoyable. In particular, I would like to acknowledge my colleagues who regularly joined me in the exhausting trenches of high-pressure research. Thank you to Matt Lucas, who taught me an innumerable number of useful things and shaped my experience at Caltech. Thank you to Jorge Munoz and Lisa Mauger, without whom the work in this thesis would not have been possible. I would like to acknowledge and thank the scientists and staff at the National Synchrotron Light Source and the Advanced Photon Source who made the experimental portions of this thesis possible. To my parents, Mike and Treo, you have been there for me from day one and I will never forget all that you have done to help me. I cannot imagine better parents than you two. To my sisters, Amber and Jenna, I could not ask for better siblings. To Grandpa Ab, thank you for your generous financial support and for your surprisingly insightful advice on practical applications for my research. To Nana and Popper and Grandma Jackie, your love and support have shaped me. Lastly, I would like to thank my wife, Jaime Marie Winterrose. Your patience, flexibility, willingness to let me work whatever hours I needed to, goodwill, and love are what got me through to this point. My life could not be complete without you. Financial support for the work presented in the pages that follow was provided by the CarnegieDOE Alliance Center, the California Institute of Technology, the John and Ursula Kanel Fellowship, the Department of Energy, and the National Science Foundation. vi Abstract The Invar effect has remained at the forefront of materials research since Charles-Édouard Guillaume discovered the vanishing thermal expansion of Fe-Ni alloys in 1897. More recently, a pressure-induced Invar effect was discovered in Fe-Ni alloys, and the relationship between classical and pressureinduced Invar phenomena has added complexity to the century-old struggle to comprehend the microscopic origins of Invar behavior. In this thesis I present our recent discovery of pressure-induced Invar behavior in Pd3 Fe with the ordered L12 structure. Nuclear forward scattering measurements show that the ferromagnetic ground state in Pd3 Fe is destabilized with pressure, collapsing around 10 GPa (V/V0 =0.96) to a lowspin magnetic state. From high-pressure synchrotron x-ray diffraction measurements we find a large volume collapse at ambient temperature to accompany the collapse of ferromagnetism. After the volume collapse there is a significant increase in the bulk modulus. Using nuclear resonant inelastic x-ray scattering to study the 57Fe phonon partial density of states (PDOS) at high pressures, we find the pressure-induced magnetic transition to cause an anomalous relative softening of the average phonon frequency. Heating our sample to 650 K in a furnace at a pressure of 7 GPa, synchrotron xray diffraction measurements reveal negligible thermal expansion from 300 to 523 K, demonstrating pressure-induced Invar behavior in Pd3 Fe. Density functional theory calculations identify a ferromagnetic ground state in Pd3 Fe with large moments at the Fe sites. These calculations show that the application of pressure counteracts the band-filling effect of Pd. By tuning the position of the top of the 3d band with respect to the Fermi level, pressure-induced Invar behavior resembles classical Invar behavior that is controlled by chemical composition. This insight marks the first step towards a unification of our understanding vii of classical and pressure-induced Invar behavior. Pressure drives the majority-spin t2g antibonding electronic states closer to the Fermi level. The transition to the low-spin state occurs as these t2g states move across the Fermi level, transferring charge to the minority-spin eg nonbonding electronic states. This charge transfer reduces the internal electronic pressure in the material, giving a volume reduction in the low-spin state. The movement of the t2g states with increasing pressure results in a greater number of states at the Fermi level, increasing screening efficiency and softening the first nearest-neighbor Fe-Pd longitudinal force constants in the low-spin state. The measured and calculated magnetic transition pressures differ significantly, despite sharing similar elastic properties in both the ferromagnetic and low-spin states. The magnitude of the disagreement between theoretical and experimental magnetic transition pressures suggests a spin-disordered state exists at high pressures in Pd3 Fe. A shape discrepancy between the calculated and measured high-pressure Fe PDOS suggests significant short-range spin correlations exist in this spin-disordered state. viii Contents Acknowledgements iv Abstract vi 1 Introduction 1 1.1 Invar History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.1 Weiss 2-γ Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.2 Noncollinear Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.1.3 Disordered Local Moment Model . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.1.4 t2g -eg Charge Transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.1.5 Pressure-Induced Invar Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2 Experimental Foundations 16 2.1 Pressure Generation: The Diamond Anvil Cell . . . . . . . . . . . . . . . . . . . . . 16 2.2 Equation of State Characterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.1 Weiss-Like Equation of State . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.2 Synchrotron X-Ray Diffraction . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Nuclear Resonant Scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.3.1 Nuclear Forward Scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.2 Nuclear Resonant Inelastic X-Ray Scattering . . . . . . . . . . . . . . . . . . 33 2.3 3 Theoretical Foundations 3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 35 ix 3.2 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2.1 Born-Oppenheimer Approximation . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2.2 Hohenberg-Kohn Theorems: Charge Density as Fundamental Quantum Variable 39 3.2.3 Kohn-Sham Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.2.3.1 Kohn-Sham Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . 42 3.2.3.2 Solving the Kohn-Sham Equations Iteratively . . . . . . . . . . . . . 43 Exchange-Correlation Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2.4.1 Local Density Approximation . . . . . . . . . . . . . . . . . . . . . . 48 3.2.4.2 Generalized Gradient Approximation . . . . . . . . . . . . . . . . . 49 3.2.4 3.2.5 . . . . . . . . . . . . . . . . . . . . 50 3.3 First-Principles Lattice Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.4 Quantum Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.4.1 Variational Quantum Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . 59 3.4.2 Diffusion Quantum Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.4.3 Quantum Monte Carlo Example: C2H4 . . . . . . . . . . . . . . . . . . . . . 61 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.5 The Projector Augmented Wave Method 4 Pressure-Induced Invar Behavior in Pd3Fe 65 4.1 Experimental Methods and Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2 Density Functional Theory Calculations . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5 Iron atom dynamics across the pressure-induced Invar transition in Pd3Fe 80 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.2.1 Nuclear Forward Scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.2.2 Nuclear Resonant Inelastic X-Ray Scattering . . . . . . . . . . . . . . . . . . 84 x 5.3 5.4 5.2.3 Born-von Kármán Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.2.4 Density Functional Theory Calculations . . . . . . . . . . . . . . . . . . . . . 86 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.3.1 Average Phonon Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.3.2 Force Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.3.3 Electronic Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.3.4 Lattice Dynamics Across the Invar Transition . . . . . . . . . . . . . . . . . . 90 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6 Concluding Remarks 96 A Flux, Brightness, and Brilliance 102 B Ab Initio Temperature Dependent Enthalpy of Reaction in Li-Mg-Si System 104 Bibliography 113 xi List of Figures 1.1 Thermal expansion coefficient in Fe-Pd system (adapted from [2]). . . . . . . . . . . . 2 1.2 Density functional theory energy-volume curves in Invar alloy. . . . . . . . . . . . . . 5 1.3 Noncollinear magnetic ground state in Fe-Ni (reproduced from [22]). . . . . . . . . . . 6 1.4 Noncollinear energy-volume curve (adapted from [22]). . . . . . . . . . . . . . . . . . . 7 1.5 Spontaneous volume magnetostriction and magnetic moment magnitude in Fe-Pd (reproduced from [32]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.6 Electronic density of states for Fe-Pt alloys (reproduced from [32]). . . . . . . . . . . . 10 1.7 Geometry of t2g and eg orbitals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.8 t2g -eg scenario for Invar behavior. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1 Merrill-Bassett diamond anvil cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Details of Merrill-Bassett diamond anvil cell. . . . . . . . . . . . . . . . . . . . . . . . 18 2.3 Diamond anvil cell mounted on Advanced Photon Source beamline. . . . . . . . . . . 19 2.4 Density functional theory results demonstrating the Weiss 2-γ model. . . . . . . . . . 20 2.5 Increase in x-ray brilliance during the twentieth century. . . . . . . . . . . . . . . . . . 24 2.6 The Advanced Photon Source at the Argonne National Laboratory. . . . . . . . . . . 25 2.7 High-pressure x-ray diffraction data from the National Synchrotron Light Source at Brookhaven National Laboratory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.8 X17C beamline at the National Synchrotron Light Source. . . . . . . . . . . . . . . . 27 2.9 57 Fe nuclear energy levels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.10 Nuclear forward scattering spectra shape when no hyperfine field is present. . . . . . . 32 xii 2.11 NFS time spectra resulting from the collective, coherent excitation of two nuclear energy levels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.1 Logic of density functional theory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2 Iteration to self-consistency to solve the Kohn-Sham equations. . . . . . . . . . . . . . 45 3.3 The local density approximation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.4 The projector augmented wave method. . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.5 Potential energy of placing a hydrogen atom at various distances above the surface of a C2 H4 molecule. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.6 Electron density contour of a hydrogen atom above an ethylene molecule. . . . . . . . 63 4.1 L12 crystal structure for Pd3 Fe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2 Volume-pressure data obtained from synchrotron x-ray diffraction measurements. . . . 67 4.3 Nuclear forward scattering spectra from Pd57 3 Fe. . . . . . . . . . . . . . . . . . . . . . 69 4.4 Diamond anvil cell furnace on beamline. . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.5 Pressure-induced Invar effect in Pd3 Fe. . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.6 Spin directions at Fe sites in the calculated magnetically ordered states in Pd3 Fe. . . 72 4.7 Energy-volume data from DFT LSDA calculations on Pd3 Fe. . . . . . . . . . . . . . . 73 4.8 Charge at the Fe site from DFT calculations. . . . . . . . . . . . . . . . . . . . . . . . 74 4.9 Magnetic moments at the Fe site from DFT calculations. . . . . . . . . . . . . . . . . 75 4.10 Experimental pressure-volume data together with DFT-calculated equation of state. . 76 4.11 Electronic density of states for alloys from the Fe-Pd system. . . . . . . . . . . . . . . 77 4.12 Electronic density of states in Pd3Fe at various volume compressions. . . . . . . . . . 78 5.1 Low-temperature nuclear forward scattering spectra from Pd3 57Fe . . . . . . . . . . . 83 5.2 Nuclear resonant inelastic x-ray scattering measurements on L12-ordered Pd3 Fe. . . . 85 5.3 Measured and calculated average vibrational frequencies in Pd3 Fe. . . . . . . . . . . . 88 5.4 Fe force constants from fit to a Born-von Kármán model. . . . . . . . . . . . . . . . . 89 5.5 Electronic density of states at the Fe site in Pd3 Fe. . . . . . . . . . . . . . . . . . . . . 91 xiii 5.6 Fe partial phonon density of states from NRIXS measurements and DFT calculations. 92 5.7 Phonon dispersion curves for Pd3 Fe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.8 Atomic displacements associated with the high-frequency mode at the X-point in Pd3 Fe. 94 6.1 t2g and eg orbitals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 Crystal structures for Li2MgSi calculation. . . . . . . . . . . . . . . . . . . . . . . . . 108 B.2 Calculated phonon density of states for Li2MgSi calculation. . . . . . . . . . . . . . . 109 B.3 Energy-separation data generated for the H atoms in the H2 molecule. . . . . . . . . . 110 B.4 Temperature dependent enthalpy of reaction for the reaction. . . . . . . . . . . . . . . 112 98 xiv List of Tables 4.1 Parameters obtained from DFT calculations of energy versus volume, fitted to the Murnaghan equation of state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Force constants projected along bond directions (longitudinal) and average perpendicular (transverse) at V/V0 =1.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 72 86 Parameters obtained from DFT (GGA-PBE)calculations on Li2 MgSi. . . . . . . . . . 107 1 Chapter 1 Introduction A new era in the century-old quest to understand Invar behavior began in 2001 with the discovery of the pressure-induced Invar effect in Fe-Ni alloys at compositions away from the classical Invar composition [1]. For the first time, it was shown that pressure could transform a material that expanded normally at ambient pressure into an Invar material. In this thesis, the first steps are taken beyond this initial discovery, showing pressure-induced Invar behavior to exist outside of the Fe-Ni system and studying the phenomenon for the first time using first-principles calculations and nuclear resonant scattering techniques. Amidst an explosion of experimental Invar-related discoveries that began in the 1950’s, sparked by a rapidly growing understanding of itinerant electron magnetism, Kussmann and Jessen [2] discovered Invar behavior in Fe70Pd30 in 1962. Figure 1.1 shows the thermal expansion coefficient (α) as a function of Pd atomic percent in the Fe-Pd system published by Kussmann and Jessen. At 30at.%Pd the thermal expansion coefficient goes rapidly to zero. This vanishing thermal expansion at a narrowly defined chemical composition is the classical Invar effect. This thesis will chronicle the experimental and theoretical work behind the discovery of pressure-induced Invar behavior in Pd3 Fe, a material far from the classical Invar composition, as shown in Fig. 1.1. A necessary precursor to the discovery of pressure-induced Invar behavior was the development of mature diamond anvil cell technology and powerful synchrotron x-ray radiation sources. Chapter 2 will introduce these developments, along with two powerful scattering techniques based on a peculiar property of a rare isotope of Iron (57 Fe). 2 Figure 1.1: Thermal expansion coefficient as a function of Pd content, reproduced from [2] with modifications The development of accurate and efficient first-principles computational methods had profound consequences for the whole field of Invar studies. Today, results from these methods lead the way for researchers striving to pin down the ever-elusive origins of classical Invar behavior. With the arrival of pressure-induced Invar studies, ab initio techniques have been pushed even further to the fore. High-pressure techniques are, after all, among the most powerful means to probe the electronic structure of materials [3]. Changing the unit cell volume in a material typically has dramatic consequences for electronic behavior, often fundamentally altering the nature of a material. The techniques introduced in Chapter 3 will allow for the direct interrogation of the electronic structure changes associated with pressure-induced Invar behavior that will be presented in Chapters 4 and 5. Iron is an amazing element, and a brief bird’s eye look at its significance is warranted before we zoom-in on its Invar-relevant details. Iron’s unique and surprising properties shape our cosmos, our 3 planet, our technology, and our bodies in ways both subtle and profound. Forged in the heart of stars, the position of iron near the maximum of the nuclear binding energy curve means that the its creation marks the culmination of the life of a star. Unable to extract energy from further fusion reactions, the nuclear fires go out and the star’s violent death throes ensue. The molecular cloud from which our planet formed was enriched by the iron end product of two earlier generations of stars, making iron the main constituent of the Earth by mass. The creation of tools and weapons made from iron heralded the rise of modern governmental and societal forms from earlier tribal societies. Our civilization continues to rely on an iron technological backbone as we begin to reach for the stars. With its prevalence and centrality likely to continue into the indefinite future, an ever-deepening understanding of the properties of iron and its alloys is of central importance to the continued progress of humanity. This thesis aims to be a further step in the long quest to understand and exploit the properties of iron. 1.1 Invar History The story of the Invar effect begins in the 1890’s with Charles-Édouard Guillaume, at the International Bureau of Weights and Measures in France, searching for a viable means of extending the standard for length measurement developed in Europe to the rest of the world. At the time, a meter was officially defined to be the length of a platinum-iridium bar kept in Europe. Materials typically expand when heated due to anharmonicities in the atomic potential wells. Platinum-iridium has particularly low thermal expansion and so was an ideal length standard. Unfortunately, platinumiridium was too costly to allow for broad replication; a less expensive replacement was needed. In his quest, Guillaume carried out a systematic investigation of the Fe-Ni family of alloys. In 1897 he happened upon the Fe64Ni36 alloy, which he discovered did not expand when heated over a large temperature range around room temperature [4]. This was exactly what Guillaume, and the world, was looking for; an inexpensive alloy showing high dimensional stability. He named the material Invar because of its invariant volume under heating. Guillaume was awarded the 1920 Nobel Prize in Physics for his discovery, beating out Albert Einstein, the favorite to win that year. Guillaume’s 4 1920 Nobel Prize for the discovery of the Invar effect in Fe64Ni36 was the first, and remains the only, Nobel Prize awarded for an achievement in metallurgy. Guillaume’s discovery was immediately put to use in numerous precision devices, and Invar behavior continues to be exploited today in an ever greater variety of ways. It has found uses in everything from toasters to shadow masks in computer screens to satellites, and in many other items which we do not have space to list here. Since Guillaume’s initial work, Invar behavior has been discovered in the Fe-Pt, Fe-Pd, Fe-Mn, Ni-Mn, Co-Mn, Co-Fe, Co-Ni, Cr-Fe, Cr-Mn, Fe-Ni-Cr, Fe-Ni-Mn, Fe-Ni-Co, Ni-Co-Mn, and Co-FeCr families [5], to name a few. In general, Invar materials involve at least one 3d transition metal element [5], and exhibit thermal expansion anomalies. A wide variety of other anomalous physical behaviors have also been linked to the Invar effect in various materials. Today “Invar behavior” has come to mean a broad range of anomalies in magnetic and mechanical properties. Some of these are anomalies in molar volume, elastic moduli, electrical resistivity, heat capacity, magnetization and high field susceptibility, and others [5]. 1.1.1 Weiss 2-γ Model The theoretical framework which has dominated Invar studies for the last several decades was put forward by R. J. Weiss in 1963 [6]. In the Weiss 2-γ (or 2-state) model of Invar behavior, two nearly degenerate magnetic states are postulated to exist in Invar materials. A state (denoted γ2 ) with large magnetic moments and high volume is the Invar ground state. Slightly higher in energy is a state (denoted γ1 ) with smaller magnetic moments (ordered antiferromagnetically) and smaller ground state volume. At the Invar composition, the low-spin, low-volume state is thermally excited at temperatures around room temperature, canceling the normal lattice expansion on heating caused by the lattice vibrations. The Weiss 2-γ model gained in popularity in the 1990’s when density functional theory (see Chapter 2) calculations revealed the presence of a low-spin (LS) magnetic state nearly degenerate with the ferromagnetic ground state in Fe-Ni [7, 8, 9], Fe-Pt [10, 11], and Fe-Pd [12, 13] Invar alloys. The LS state, typically lying on the order of 1 mRy (150 K) above the ferromagnetic (FM) ground 5 state, was identified with Weiss’s low-volume antiferromagnetic state. Figure 1.2 shows the typical Invar situation revealed by these ab initio calculations, with a large-volume ferromagnetic state lying just below a low-volume, low-spin state. Temperature or pressure depopulates the ferromagnetic ground state in favor of the low-spin state, causing the Invar anomaly. Figure 1.2: Density functional theory calculations reveal a smaller volume low-spin state existing at a slightly higher energy than the ground state high-volume ferromagnetic state near the classical Invar composition in Fe-Pt. Experimental support for the Weiss 2-γ model has come from high-pressure experiments showing a magnetic collapse to occur in classical Invar alloys at high pressures [14, 15, 16], as well as from experimental studies showing sudden stiffening of the bulk modulus at high pressures [17, 18, 19, 20]. However, the magnetic inhomogeneities predicted by the Weiss 2-γ model have never been detected experimentally [21]. This, along with other concerns [22], led to the development of a new set of ideas based around noncollinear magnetic states hypothesized to exist in Invar alloys. 6 1.1.2 Noncollinear Model Wang, et al. [23], were the first to calculate a noncollinear ground state in Fe65Ni35 Invar. They found a noncollinear spin alignment, in which Fe magnetic moments relaxed into canted orientations relative to the average magnetization in the material (see Fig. 1.3), to be the preferred magnetic state when Fe atoms were surrounded by less than three Ni nearest neighbors [23]. Schilfgaarde, et al. [22], formulated a theory of Invar behavior based on the noncollinear ground state in Fe65 Ni35. They found that allowing for noncollinearity resulted in an anomaly in the Grüneisen constant, leading to the thermal expansion anomaly characteristic of Invar behavior in Fe65 Ni35. Noncollinear states have also been stabilized in Fe-Pt Invar alloys [24]. Figure 1.3: Noncollinear magnetic ground state in an Fe-Ni supercell at an atomic volume of 71.9 a.u. reproduced from [22]. The calculations were carried out in the local spin-density approximation using the relativistic linear muffin-tin orbitals method. Red arrows indicate Fe moments, blue arrows indicate Ni moments. Relative arrow magnitude gives relative local magnetic moment amplitudes. The average Fe(Pd) moment at this volume was approximately 1.5(0.5) µB . Relaxing the ferromagnetic constraint in the calculations smeared out the sharp high-spin and low-spin curves, leading to a continuous transition from a noncollinear state of mostly ferromagnetic character in the ground state to an increasingly noncollinear magnetic state at smaller volumes, as shown in Fig. 1.4. These noncollinear calculations gave an improved prediction for the volume 7 change associated with the magnetic transition under pressure and explained elastic anomalies that were difficult for the Weiss 2-γ model to handle [25]. Figure 1.4: Energy-volume curves reproduced from [22] (with modification) for various magnetic states in Fe-Ni. The calculation was carried out in a 32 atom supercell containing 21 Fe atoms and 11 Ni atoms randomly distributed on an fcc lattice. However, neutron scattering with polarization analysis failed to identify noncollinear magnetic states in Fe65Ni35 [26]. Instead the measurements showed Fe65Ni35 to be a conventional ferromagnet in its ground state. The incorrect ab initio prediction of a noncollinear ground state in classical Invar alloys is now believed to stem from the overbinding inherent in the local-spin-density approximation [27]. The reduced ground state volume in the LSDA calculations stabilizes the noncollinear states. Some high-pressure experimental studies have supported the noncollinear model [28, 29, 30], and the controversy continues to the present. 1.1.3 Disordered Local Moment Model With both the Weiss 2-γ and noncollinear models apparently failing to account for all aspects of Invar behavior, an attempt has recently been made to place Invar’s origin in the realm of typical magnetic material behavior [21, 31, 32, 33]. In this view, the separation between Invar and non-Invar materials is of a quantitative, and not a qualitative, nature. 8 A key measure of Invar behavior is the spontaneous volume magnetostriction, ωs = V0 (F M ) − V0 (P M ) , V0(P M ) (1.1) defined as the relative volume change between the FM and paramagnetic (PM) states. Most magnetic materials show some spontaneous volume magnetostriction at the disordering temperature. Invar alloys are distinguished by the large size of their ωs (10−2 is typical) [5]. Large, positive values for ωs in Invar alloys cancel the usual thermal expansion, resulting in zero thermal expansion (the Invar effect). To explain why Invar materials exhibit a particularly large ωs , Khmelevskyi, et al. [21], calculated ωs for a series of Invar and non-Invar materials. Traditionally, ωs is poorly predicted by density functional theory methods. The DFT representation of the low-spin state as a nonmagnetic state in which magnetic moments completely vanish at the atom sites is the cause of the inaccuracy. Local, fluctuating moments are believed to exist in many materials above their magnetic ordering temperatures, and this has consequences for the spontaneous volume magnetostriction. Khmelevskyi was able to obtain accurate predictions for ωs using the disordered local moment (DLM) method [34] to account for fluctuating local moments in the paramagnetic state. In the DLM method a binary alloy is treated as a pseudo ternary alloy by treating the spin-up Fe magnetic moments separately from the spin-down Fe magnetic moments atoms [31]. In the ground state, all spins are aligned. At the Curie temperature, half of the Fe moments point in one direction and the other half of the Fe moments point in the opposite direction, giving the paramagnetic state. Figure 1.5a shows the calculated spontaneous volume magnetostriction for chemically disordered Fe-Pd alloys from [32]. Fe concentration increases towards lower valence electron concentration (towards the left) in Fig. 1.5. ωs is at a maximum at the Invar composition (Fe70 Pd30), and drops rapidly with increasing Pd content. The calculated ωs overestimates the experimental values by approximately 20% [5, 32]. Maximum spontaneous volume magnetostriction coincides with maximum change in local Fe magnetic moment magnitude, shown in Fig. 1.5b. This is understood by recalling that spin polarization increases valence electron kinetic energy. This kinetic energy increase 9 is counterbalanced by an increased volume in the spin-polarized state. In this manner, larger local moments correspond to larger ground state volumes [32]. According to the view of the DLM model, the large change in local Fe magnetic moment magnitude when going from the FM to the DLM state is the underlying cause of the large ωs at the Invar composition. Figure 1.5: (a) Spontaneous volume magnetostriction (ωs) and (b) local magnetic moments as a function of valence electron concentration at the Fe site for the Fe-Pd system, reproduced from [32] with modification. The electronic structure situation is illustrated in Fig. 1.6. The key finding in these DLM-Invar studies was that the change in the local magnetic moment magnitude became maximal when the top of the majority 3d band was just below the Fermi level (Ef ) in the ferromagnetic state [32]. Magnetic 10 disorder in the DLM state then caused changes (of dynamical origin) in the electronic density of states, with the highest energy majority spin states moving across the Fermi level, causing a maximal change in local moment magnitude. With decreasing Fe concentration the top of the majority 3d band moved away from the Fermi level (Fig. 1.6b), stabilizing the strong ferromagnetic state and leading to a reduced moment change on spin disordering, explaining the disappearance of Invar behavior away from a narrowly defined composition range. Figure 1.6: Calculated spin-polarized electronic density of states in (a) Fe70 Pt30 Invar alloy and (b) Fe50Pt50 alloys, reproduced from [32]. This model has had some success with the Fe-Pt and Fe-Pd classical Invar alloys, and claims have been made that classical Invar behavior is now understood in these alloys [21]. The DLM model gave controversial results for Fe-Ni, however, leading to the formulation of a theory of Invar 11 behavior in Fe65Ni35 based on magnetochemical effects [35, 27]. In this model the Invar anomalies arise due to local chemical order in disordered Fe-Ni alloys. Another weakness of the DLM picture is that it does not account for short range magnetic correlations in the spin-disordered state [36], which are known to play a role in Invar materials [37]. 1.1.4 t2g -eg Charge Transfer Electronic structure calculations have provided a detailed picture of the charge transfers around the Fermi level in Invar materials. Figure 1.7 shows the majority-spin electronic densities of states (DOS) at the Fe site for L12 Pd3 Fe, decomposed into states of t2g and eg character. In an fcc lattice containing transition metal alloys, the t2g orbitals form strong ddσ bonds with first-nearest neighbor (1NN) atoms, owing to their large charge density in the [110] direction. This causes the t2g DOS to split into a high-energy sub-band with anti-bonding character and a low-energy sub-band with primarily bonding character. The eg orbitals, with large charge densities between the more distant 2NN atoms along the [100] direction, form weaker ddπ bonds. They are energetically intermediate between the bonding and anti-bonding t2g states, as shown in Fig. 1.7. In Invar alloys in the ferromagnetic state, nonbonding minority-spin eg states just above the Fermi energy are nearly degenerate with t2g majority-spin antibonding states just below the Fermi energy. Kaspar and Salahub [38] were the first to recognize this fact, based on molecular-orbital cluster calculations on Fe-Ni metal clusters, and use it to explain Invar behavior. Entel, et al. [8], generalized Kaspar and Salahub’s earlier results to L12 ordered Fe3 Ni. A simplified schematic of the situation is given in Fig. 1.8a. The basic idea was that as charge is transferred (with pressure or temperature) from the t2g majority-spin antibonding states to the eg minority-spin nonbonding states (green arrow in Fig. 1.8), the internal electronic pressure in the material is reduced, giving a volume reduction that cancels the usual thermal expansion with temperature caused by atom vibrations in anharmonic potential wells. This t2g -eg scenario has obvious resonance with the Weiss 2-γ model (Section 1.1.1), and the t2g eg scenario has been cited repeatedly as support for the Weiss 2-γ model in the literature [20, 36, 39]. 12 Figure 1.7: (upper) Geometry of the t2g (left) and eg (right) orbitals on an L12 lattice. Blue marks regions of large charge density. (lower) Majority-spin electronic density of states at the Fe site in Pd3 Fe decomposed into states of t2g and eg character. Ef =0 eV. Weiss’s high-spin, high-volume state is associated with the occupation of t2g antibonding orbitals, while the transition to Weiss’s low-spin, low volume state is associated with the transfer of charge from the t2g antibonding states to the eg nonbonding states. Nevertheless, the importance of the relative occupation of the t2g and eg orbitals has significance beyond the Weiss 2-γ picture. Calculations allowing noncollinear spin alignment (Section 1.1.2) have revealed many nearly degenerate collinear and noncollinear magnetic states around the ground state volume in Invar alloys [40]. The origin of these degenerate states has been traced to the near degeneracy of the t2g antibonding and eg nonbonding states in Invar alloys [41, 42, 43]. Changes in temperature or pressure cause charge transfer between the t2g antibonding and eg nonbonding orbitals near the Fermi level, causing transitions between nearly degenerate collinear and noncollinear states [43]. 13 Figure 1.8: Essential features of the electronic density of states involved in the t2g -eg scenario for the (a) high-spin, and (b) low-spin magnetic states in Invar materials. The green arrow indicates the charge transfer responsible for Invar behavior with increasing temperature or pressure. The t2g -eg picture adds additional microscopic insight to the noncollinear Invar model. The position of the t2g states relative to the Fermi level is also of central importance for the theory of Invar behavior based on the disordered local moment method [32] (Section 1.1.3). A decrease in Fe content away from the Invar composition causes the majority-spin antibonding t2g levels to move away from the Fermi level. This leads to a smaller change in local magnetic moment magnitude upon spin disordering, and so a smaller spontaneous volume magnetostriction and no Invar anomaly at ambient pressure in these strong ferromagnets. Once again, it is the relative position of the antibonding t2g majority-spin states that determines the presence or absence of Invar behavior. A full t2g -eg analysis of the DLM results has yet to be carried out. 14 In addition, the t2g and eg scenario developed for Invar behavior has provided insight into a wide range of distinctive Fe alloy behavior beyond the Invar anomaly. These include anti-Invar behavior, the martensitic transformation, and phonon anomalies [43]. Polarized neutron diffraction measurements failed to find the large changes in t2g and eg electronic occupations with temperature predicted in Fe65 Ni35 [44] and Fe72Pt28 [45, 46]. These measurements and their implications are, however, not without controversy [35, 36]. High temperature paramagnetic neutron scattering on Fe-Ni alloys do appear to support the t2g -eg scenario for the origin of the Invar effect [36, 37]. The t2g -eg scenario provides a powerful theoretical framework within which Invar behavior, and many other characteristic properties of Fe and its alloys, can be understood. The explanatory structure the t2g -eg picture provides has proven useful in diverse Invar materials, and allows microscopic insight into the leading theories of Invar behavior. A clear experimental picture with regards to the t2g -eg scenario awaits additional measurements in the future. 1.1.5 Pressure-Induced Invar Effect To this point, pressure-induced Invar behavior has shown the same complexity of behavior that has characterized classical Invar alloys. The initial discovery [1] was interpreted in terms of the noncollinear model. The work presented in this thesis spans the theoretical landscape, incorporating aspects of the Weiss 2-γ model, noncollinear, disordered local moment, and the t2g -eg scenario. Studying pressure-induced Invar behavior is a uniquely powerful means to come to grips with classical Invar behavior because, with the pressure-induced case, we are able to start with a non-Invar metal and transform it slowly into an Invar metal through the application of pressure. By identifying the changes pressure initiates that transforms the material into an Invar alloy, we can come to understand the essential microscopic origins of Invar behavior that have eluded researchers for over a century. After the discovery of classical Invar behavior in 1897 it took many decades before the burst of discoveries identifying Invar behavior in many other alloy systems took place. It is a mark of 15 the accelerated rate of progress we currently find ourselves experiencing that we are now on the cusp of a similar explosion of discoveries in the field of pressure-induced Invar studies. Our work on Pd3 Fe presented in this thesis marks the first move beyond the Fe-Ni alloys in this regard. Many other systems exhibiting pressure-induced Invar behavior await, just beyond our current view, to be discovered. May the work presented in this thesis guide future researchers embarked on this exciting journey of discovery. 16 Chapter 2 Experimental Foundations In this chapter we will explore developments at the intersection of nuclear condensed matter physics and high-pressure research that have recently come together to open a powerful new window on the solid state. 2.1 Pressure Generation: The Diamond Anvil Cell Pressure (P) is defined as, P = F , A (2.1) with F the normal force and A the area over which the force is applied. From this simple equation it is immediately apparent that to maximize pressure one should minimize A and maximize F. With its extreme hardness and incompressibility, diamond is the material of choice for this task. With diamond’s broad transparency to electromagnetic radiation (with the notable exceptions of ultraviolet and soft x-rays) pressure cells built around diamond anvils are the key tool for the highpressure researcher. In simplest terms, a diamond anvil cell (DAC) consists of two diamonds and the means to push them together. Many varieties of DAC are in use today. Figure 2.1 shows a variation of the MerrillBassett type diamond anvil cell [47, 48] in its assembled (top left) and disassembled (bottom) form. The advantage of this so called ‘Tel-Aviv’ type diamond anvil cell is its small size (note the stack 17 Figure 2.1: Merrill-Bassett diamond anvil cell used in this thesis. Upper: fully assembled DAC and stack of U.S. quarters for size comparison. Lower: disassembled DAC. (Photo by Chen Li.) of quarters in Fig. 2.1) and relative simplicity of design and use. This simple design, consisting of two opposed, identical backing plates, is also the source of the major drawback with this DAC. At moderately high pressures the backing plates begin to deform, making higher pressures unattainable. Figure 2.2 shows the essentials of how the Tel-Aviv type diamond anvil cell operates. The diamonds are affixed to two backing plates, with culets facing each other. A thin piece of metal, known as a gasket, rests between the two diamond culets. A small hole (∼200 µm) is drilled in the center of the gasket. The sample is placed in the center of this hole. The sample chamber, formed around the sample by the diamond culets and drilled gasket, facilitates the use of a pressure transmitting medium (typically silicone oil), which surrounds the sample and ensures hydrostatic pressure is transmitted to the sample. The use of a metal gasket also helps to prevent diamond fracture at high pressures. Stainless steel is often the simplest choice for gasket material. In measurements that are particularly sensitive to Fe content, Rhenium can be used in its place. If equatorial access is a requirement for the measurement of emitted photons, a beryllium gasket is used. The diamonds are forced together by the tightening of six screws. The sample, resting between the diamond culets in the center of the gasket hole, is compressed as the diamonds are drawn together. The pressure within the sample chamber must be known to high accuracy. The vital breakthrough in this regard came with the realization that the fluorescence lines of ruby (Al2 03:Cr) shift linearly with pressure [49]. This insight was arrived at by measuring the position of the ruby fluorescence 18 Figure 2.2: Cutaway of a Merrill-Bassett, Tel-Aviv type diamond anvil cell (left) with zoom-in on the critical pressure-generating region (right). (Cutaway on left created by Lisa Mauger.) lines together with x-ray diffraction peaks from NaCl. Fitting the NaCl results to an equation of state, the pressure dependence of the ruby fluorescence lines was found. In practice, a small amount of ruby is placed in the sample chamber together with the sample. The DAC is sealed and pressure generated by tightening the DAC screws. Laser light enters through the axial hole in the DAC and excites the ruby. The wavelength of the fluorescence gives the pressure in the sample chamber. Multiple measurements from different parts of the sample chamber are often taken to ensure accurate pressure calibration and to provide error bars from pressure variations. With the pressure calibrated, the DAC is mounted on the synchrotron beamline (Fig. 2.3). Synchrotron x-rays enter the DAC through the axial hole and scatter from the sample to give the desired measurement. Much of the real-time effort on the beamline during a high-pressure synchrotron experiment consists of scanning the finely focused (20 µm (horizontal) by 50 µm (vertical) at the 16ID-D APS beamline, for example) synchrotron x-ray beam in search of the gasket hole and sample. To obtain a reliable measurement, one must ensure that the beam is hitting the sample and 19 Figure 2.3: Diamond anvil cell mounted on the 16ID-D beamline at the Advanced Photon Source. (Photo by Matt Lucas.) not the gasket edge or DAC. This searching often requires a high-degree of persistence on the part of the researchers, adding pressure in the limited time for their experiments. To avoid diamond fracture and gasket failure at high pressures, the two diamond culets must be positioned horizontally to one another to a high degree of accuracy. Achieving this, while only tightening one screw at a time, is an additional difficult task. It is here that much of the art in DAC work enters, with expertise often taking years to develop. 2.2 Equation of State Characterization An equation of state is a relation among the state variables that characterize the thermodynamic state of a system. In experimental high-pressure work the important state variables are typically pressure (P), specific volume (V), and temperature (T). When working with solids, temperature is often neglected in the equation of state due to the relative insensitivity of the volume of solids to temperature [50]. We are then left with relating the pressure exerted on a material to its specific volume. This P-V relationship is an isothermal equation of state. A number of isothermal equations of state are in widespread use in high-pressure research today. As we will see in Chapter 4, Pd3 Fe undergoes a large volume collapse between 10 and 15 GPa, and a subsequent decrease in compressibility above 15 GPa. Standard isothermal equations of state cannot 20 account for these dramatic high-pressure behaviors. Recently, the “Weiss-like” equation of state [20] was developed to model such magnetovolume effects, a topic to which we turn now. 2.2.1 Weiss-Like Equation of State The ideas behind the Weiss-like equation of state originate in the Weiss 2-γ model of Invar behavior (outlined in Chapter 1). In this model, proposed by Weiss in 1963 [6], two nearly degenerate magnetic states are postulated to exist in Invar materials (see Chapter 1 for an explanation of Invar). A state with large magnetic moments and high volume is the Invar ground state. Slightly higher in energy is a state with vanishing magnetic moments and smaller ground state volume. This is shown in Fig. 2.4a. At the Invar composition, the low-spin, low-volume state is thermally excited at temperatures around room temperature, canceling the normal lattice expansion on heating caused by the lattice vibrations. Figure 2.4: (a) Density functional theory results showing a low-volume paramagnetic (low-spin) state existing at a slightly higher energy than the ground state high-volume ferromagnetic state in a typical Invar alloy. These two states, first proposed to exist in Fe-Ni by Weiss [6], form the backbone of the 2-γ model of Invar behavior. (b) As Fe content decreases, the ferromagnetic state is stabilized and Invar behavior disappears. While iron-rich fcc based Invar alloys are characterized by the existence of nearly degenerate states with differing magnetic moments and volumes (Fig. 2.4a), decreasing Fe concentration stabi- 21 lizes the ferromagnetic state (Fig. 2.4b), and suppresses classical Invar behavior. This means that away from the Invar composition, the high-spin, ferromagnetic state is fully occupied. The Weisslike equation of state models this low-pressure, high-spin region with a conventional Murnaghan [51] equation of state,  −(1/B00 ) 0 P V (P ) = V0 1 + B0 , B0 (2.2) 0 with B0 the zero-pressure isothermal bulk modulus, B0 the pressure derivative of the bulk modulus, and V0 the ground state volume. This can be derived by expanding the bulk modulus to first order 0 in P about P=0 (making the assumption that B0 is independent of pressure), 0 B = B0 + B0 P, (2.3) and recalling that, B = −V  ∂P ∂V  . (2.4) T Inserting Eq. 2.4 into Eq. 2.3 and integrating gives Eq. 2.2. As pressure increases, the ferromagnetic state becomes less favorable and electrons are eventually driven to the low-spin state. After the electrons have been fully driven to the low-spin state, the thermal expansion is again normal. The Weiss-like equation of state treats this high-pressure region with a second Murnaghan equation of state. The region between these two normal compression regions is more complicated. In this region, the electrons are in the process of transition to the low-spin state. This transition region is marked by a rapid volume collapse, which is modeled as a weighted-sum of the high-pressure and low-pressure Murnaghan equations of state. The weighting factor is given by a thermodynamic occupancy factor, n(P ). To find the form of the occupancy factor, we recall that the population of high-spin electronic states is given by the Fermi-Dirac distribution, 22 n(E) = 1 , 1 + e−∆E/(kB T ) (2.5) with, ∆E = EHS − ELS . (2.6) Expanding ∆E to first order in P about P = 0, and making the assumption that the electronic states have energies that depend linearly on pressure, we obtain ∆E = ∆E0 − ∂∆E P. ∂P (2.7) We substitute Eq. 2.7 into Eq. 2.5, rearrange the terms, defining, PT = ∆E0 ∂P , ∂∆E (2.8) ∆P = kB T ∂P , ∂∆E (2.9) to obtain, n(P ) = 1 1 + e(P −PT )/∆P . (2.10) Here PT is the pressure of the transition between the high-spin (HS) and low-spin (LS) regions, and ∆P is the pressure ranger over which the magnetic transition takes place. Putting all of these pieces together, we arrive at the Weiss-like equation of state, 23 V (P ) =  1 1 + e(P −PT )/∆P   + 1−    0 V0,HS 1 + B0,HS 1 1 + e(P −PT )/∆P  P B0,HS   0 V0,LS 1 + B0,LS 0 −(1/B0,HS ) P B0,LS   0 −(1/B0,LS )  . (2.11) We will see the power of Eq. 2.11 to model pressure-induced rapid volume collapse regions in Chapter 4. To put our equation of state to use, we need to determine the pressure and volume of a material while under compression. As mentioned above, the fluorescence of ruby technique [49] gives the pressure of the sample in the diamond anvil cell. We are then left with the need to determine the volume of the sample under pressure. Modern third-generation synchrotron radiation sources allow very efficient characterization of the specific volume of a material by means of x-ray diffraction. 2.2.2 Synchrotron X-Ray Diffraction X-rays, with wavelengths on the order of typical interatomic distances in crystals, are ideal for diffraction studies. These studies reveal crystal structure and unit cell volume. The small sample sizes required by the demands of high-pressure research make synchrotron x-ray sources the tool of choice for probing the properties of materials at high pressures. Magnetic fields cause charged particles, such as electrons, to move in circular orbits. This centripetal acceleration causes the emission of electromagnetic radiation, in this case known as cyclotron radiation. Synchrotron radiation is the name given to such radiation when it is emitted by charged particles accelerated to relativistic speeds. A synchrotron radiation source is essentially a large evacuated ring around which electrons circulate. As the electrons orbit, they continuously emit synchrotron radiation. Beginning in the late 1970’s synchrotron rings began to be built exploiting specialized magnet arrangements for local control of electron acceleration, increasing brilliance by several orders of magnitude over earlier sources [52]. This strategy, using magnet arrays to increase brilliance, defined as [52], 24 Brilliance = P hotons/s mrad2 mm2 0.1%bandwidth (2.12) (photons per second within a specified energy bandwidth, normalized to the solid angle of the radiation cone and to the source area), was so successful that soon synchrotron sources were built around the use of these insertion devices. These third generation sources incorporated long, uncurved sections in which linear arrays of dipole magnets, called “undulators” were inserted. The orientation of the magnets in an undulator alternate, causing electrons to oscillate rapidly in sinusoidal paths, emitting very intense synchrotron radiation. By adjusting the spacing between the magnets the characteristics of the radiation can be finely tuned to the needs of the experimentalist. Figure 2.5 shows the rapid increase in spectral brilliance the use of insertion devices has made possible. Figure 2.5: Increase in x-ray brilliance during the twentieth century. Figure adapted from [52], data from [53]. X-rays incident on a crystal interact with the electron clouds around the atoms, the electric field from the incident x-ray accelerating the charged electrons, causing re-radiation (elastic scat- 25 Figure 2.6: The Advanced Photon Source (APS) at the Argonne National Laboratory. (Upper) Aerial view of the synchrotron ring. (Lower) Arrangement of a typical beamline at the APS. Images extracted from [54]. tering) [55]. Radiation scattered by two parallel planes of atoms (separated by a distance d) show constructive interference when the wavelength of the scattered radiation (λ) is equal to the length difference traveled by radiation scattered from the different respective planes. This is Bragg’s Law [56], summarized in Eq. 2.13, 2dsinθ = λ. (2.13) Two basic strategies are routinely employed to obtain x-ray diffraction patterns. The first is 26 angle-dispersive x-ray diffraction (ADXD), which employs a monochromatic beam of x-rays and a detector that records the scattering at different angles. The synchrotron radiation is monochromized by carefully selected Bragg diffraction from a crystal (called a monochromator). The second method, Energy Dispersive X-ray Diffraction (EDXD), uses a polychromatic x-ray beam and a detector at a fixed angle that is capable of recording scattered radiation with different wavelengths. Data obtained from such a measurement, performed at beamline X17C of the National Synchrotron Light Source at Brookhaven National Laboratory, is shown in Fig. 2.7, together with the beamline environment in Fig. 2.8. The x-ray diffraction peaks are indexed (using, for example, Rietveld analysis [57]) to give the lattice parameter at each pressure. Figure 2.7: Energy dispersive x-ray diffraction patterns measured at beamline X17C at the National Synchrotron Light Source from FePd at various pressures plotted versus scattering angle (2θ). 27 Figure 2.8: X17C beamline at the National Synchrotron Light Source set up for energy dispersive x-ray diffraction. The location of the fixed angle Ge detector and diamond anvil cell (DAC) is indicated, together with the direction from which the x-ray beam enters the cell. 2.3 57 Nuclear Resonant Scattering Fe, with a 2.1% natural abundance [52], has a first nuclear excited state only 14.41 keV above its ground state. This is miniscule compared to the 847 keV gap between the ground and first excited state of the most common natural isotope of Fe, 56Fe, for example. The small energy gap between the ground and first-excited states in 57 Fe has the surprising consequence that 57Fe can be excited into its first nuclear excited state without recoil when the 57Fe atom is bound in a solid. This is the Mössbauer effect, discovered first in 191Ir by Rudolf Mössbauer in 1958 [58]. Mössbauer was awarded the 1961 Nobel Prize in Physics for his discovery. The origin of the Mössbauer effect is found in the quantum mechanics of the crystal lattice. According to quantum mechanics, a crystal lattice can only vibrate in a discrete number of vibrational modes. For a typical lattice, the characteristic energy of these vibrations (~ωE ) is several tens of meV. The Lamb-Mössbauer factor gives the probability that emission or absorption 28 of a γ-ray by the nucleus will cause no change in the nuclear state [59]. For an Einstein solid in its ground state the Lamb-Mössbauer factor takes the simple form [59, 60]   ER , fLM = exp − ~ωE (2.14) with ER the recoil energy, ωE the Einstein frequency and ~ Planck’s constant. Elements with large fLM are known as Mössbauer isotopes. 57Fe is a Mössbauer isotope, while 56Fe is not. To understand why, we can calculate the recoil energy associated with exciting the nucleus of these two isotopes. The recoil energy due to absorption or emission of a photon by a nucleus is calculated as, ER = 1 2 p2 p2 c2 , mv = = 2 2m 2mc2 (2.15) with m the mass of the recoiling nucleus, v the velocity imparted by the photon to the recoiling nucleus, p the momentum of the recoiling nucleus, and c the speed of light. For 57 Fe we obtain, ER = (14.41 keV )2 = 1.96 meV, 2 ∗ 53.0221 GeV (2.16) well below the several tens of meV typical of a vibrational mode of the Fe crystal lattice. This results in a large probability for the 57Fe nucleus to be excited recoilessly when bound in a solid. We can compare this with the case for 56Fe, where we find, ER = (847 keV )2 = 6, 885.22 meV. 2 ∗ 52.0902GeV (2.17) Clearly, the probability for recoilless emission is not great for a bound 56Fe nucleus. Other Mössbauer isotopes include 169Tm, 153Eu, 149Sm, 181Ta, and 83Kr, among others. All Mössbauer isotopes share the property of having a low-lying first nuclear excited state. 57Fe is the Mössbauer isotope of interest in this thesis. For historical reasons, we refer to radiation of identical energy by different names depending on its source of origin. The term γ-ray refers to high-energy radiation originating from transitions 29 within the nucleus. The term x-ray refers to radiation generated by electronic processes occurring outside the nucleus. Commonly, γ-rays are considered of higher energy than the x-rays that could be produced in the lab. Today this is no longer the case. Nevertheless, the naming convention remains. Detecting γ-rays from nuclear de-excitations amidst the storm of x-rays that are invariably scattered by the electrons in a sample is akin to finding the proverbial needle in a haystack. This task is made feasible by the long lifetime of the first excited state of 57Fe. The energy uncertainty of the first excited state of the 57Fe nucleus is 4.67 neV [59]. Through the uncertainty principle this gives a lifetime of 141 ns, a virtual eternity compared to the electronic scattering which takes place on the femtosecond (10−15 s) timescale. Synchrotron x-ray sources do not deliver x-rays continuously in time. Instead, electrons circle the storage ring in bunches, delivering x-ray photons in pulses. At the Advanced Photon Source, the standard operating mode (top-up) has 24 single bunches (singlets) of electrons orbiting the storage ring simultaneously [61]. Each bunch has a current of approximately 4.25 mA, with a time spacing of 153 ns between singlets. Lasting about 40 ps, a singlet transverses the insertion device and a pulse of x-rays travels towards the sample. After monochromization (which typically happens in several steps) the x-rays are incident on the sample. In Fe NRS, the incident beam is tuned to 14.41 keV to excite the 57Fe nuclei. After the initial x-ray pulse, the detector signals are switched off to allow the intense electronic scattering to pass. At around 20 ns after the initial pulse, the counting is switched on to count the delayed γ-rays as a function of time until the next x-ray pulse is delivered. It was Stanley Ruby in 1974 [62], who was the first to suggest that Mössbauer spectroscopy might be performed using synchrotron radiation. A decade after Ruby’s suggestion the first observation of nuclear resonant scattering of synchrotron radiation was made in 1985 by Erich Gerdau’s group [63]. After a further decade of effort the first successful nuclear resonant inelastic x-ray scattering measurements were demonstrated in 1995 by Seto and Sturhahn [64, 65, 66]. Today, the routine availability of third generation synchrotron x-ray sources in many parts of the world has allowed nuclear resonant scattering to become a standard tool in the materials scientist’s tool box. Combined with diamond anvil cell technology, nuclear resonant scattering gives a unique 30 and powerful window on the behavior of materials at extremes of pressure and temperature. 2.3.1 Nuclear Forward Scattering The Pauli exclusion principle tells us that many-body fermionic wavefunctions must be antisymmetric under exchange of electrons. Two electrons with parallel spin will have antisymmetric spatial wavefunctions. Electrons with oppositely oriented spins will have symmetric spatial wavefunctions. These spatial wavefunctions entail different Coulombic energy contributions to the Hamiltonian. Coulomb energy is lower when electrons have parallel spins because they stay further apart on average than when electrons have anti-parallel spins. This is the origin of the exchange energy. This principle is extended to itinerant systems of many interacting electrons by realizing that in these systems the exchange interaction takes place between many electrons simultaneously. Opposed to the reduction in Coulombic energy from unpairing electrons is the energy penalty resulting from the unpaired electrons not residing in their lowest energy state. This penalty is reciprocally proportional to the value of the electronic density of states at the Fermi level. In magnetic materials the electronic density of states at the Fermi level must be large so that the energy penalty inherent in spin-polarization is small enough to be compensated for by the exchange interaction. At the end of the 3d series the d bands are sufficiently flat (and so the electronic density of states sufficiently large) at the Fermi level to lead to spontaneous ferromagnetic order. This is not the case with 4d electron bands. The extra node in the 4d wavefunction causes the electron density to extend further from the nucleus. This implies a stronger interaction between neighboring atoms, which broadens the bands and reduces the electronic density of states at the Fermi level. Nuclear resonant scattering in the form of nuclear forward scattering (NFS) measurements provides the ability to probe the magnetic state of materials containing Mössbauer isotopes. More specifically, in NFS measurements we obtain information on the strength of the hyperfine magnetic field at the nucleus. This field is primarily produced by the Fermi contact interaction originating from shape distortion of the spin-up s electron wavefunctions caused by the exchange interaction between spin-polarized 3d electrons and s electrons at the nucleus [59]. The details of this inter- 31 action cause the hyperfine magnetic field at the 57Fe nucleus to be in a direction opposite to the lattice magnetization [59]. The hyperfine magnetic field interacts with the nuclear spins, lifting the degeneracy of the ground and first-excited state of the 57 Fe nucleus. In the presence of a hyperfine magnetic field, the ground state is split into 2 states, while the first excited state is split into 4 states, shown in Fig. 2.9. Figure 2.9: Representation of the 57Fe nuclear energy levels without (left) and with (right) hyperfine magnetic field present. Arrows indicate transitions allowed according to the dipole selection rule (∆m=0,±1). I denotes nuclear spin. The hyperfine splitting of the nuclear energy levels is on the order of neV. The pico-second synchrotron pulse at 14.41 keV has an meV band width, and so excites all nuclear hyperfine levels simultaneously. The coherent superposition of decay products with slightly different energies due to the hyperfine splitting leads to intensity beats in time, known as quantum beats. We can understand the origin of these beats by working through a simple example with a two level system, following [67]. First consider the case when there is no hyperfine splitting of the nuclear energy levels. Incident x-rays excite the first excited states of the nuclei (Fig. 2.9(left)), and the detector (an avalanche photodiode in this case) measures an intensity (I), √ I(t) = |ψ|2 ∝ | λe(−λt/2)e(−itω/2)|2 = λe−λt , (2.18) with ψ the amplitude, t the time, λ the decay constant, and ~ω the energy (E) of the excited nuclear state. Equation 2.18 is shown graphically in Fig. 2.10. 32 Figure 2.10: NFS time spectra resulting from the excitation of a single nuclear energy level. Now consider the situation when the system can be excited to two distinct, distinguishable energy levels. Decay products from these distinguishable levels give an intensity at the detector of the form, √ √ I(t) = |ψ1|2 + |ψ2|2 ∝ | λe(−λt/2)e(−itω1 /2)|2 + | λe(−λt/2)e(−itω2/2) |2 = 2λe−λt . (2.19) Equation 2.19 is of the same form as that shown in Fig. 2.10, with twice the amplitude. We once again have a simple exponential decay in time. The situation becomes more intricate when we take into account the fact that we cannot in fact know which measured photons are coming from which particular energy level. The measured decay products are coherent, and we must take the square of the sum to account for this, √ √ I(t) = |ψ1 + ψ2|2 ∝ | λe(−λt/2)e(−itω1 /2) + λe(−λt/2)e(−itω2 /2)|2 = 2λe−λt + 2λe−λt cos((ω1 − ω2)t/2) = 4λe−λt cos2 ((ω1 − ω2)t/2) Equation 2.20 is displayed graphically in arbitrary units in Fig. 2.11. (2.20) 33 Figure 2.11: NFS time spectra resulting from the collective, coherent excitation of two nuclear energy levels. The periodic oscillations in Fig. 2.11 are the quantum beats referred to above, these originating from our simple two-level system. The time separating two neighboring minima in Fig. 2.11 is given by, ∆t=h/∆E. In this way, the energy level splitting (∆E) can be extracted from the intensity modulations. When there is more than one energy splitting, as is the case when hyperfine magnetic fields are present, the quantum beats can become quite complex and advanced methods [68] and significant effort must be deployed to extract the energy splittings. When one is simply interested in knowing whether or not magnetic order exists in a material, the situation is considerably simplified. The presence of quantum beats in the NFS spectra indicates magnetic ordering, while the absence of quantum beats indicates that long-range magnetic order is absent in the material. 2.3.2 Nuclear Resonant Inelastic X-Ray Scattering Nuclear resonant scattering in the form of nuclear resonant inelastic x-ray scattering (NRIXS) provides the unique capability of measuring phonon densities of state at high pressures. In simplest terms, in NRIXS we detune the incident photons from the nuclear resonance energy in discrete steps, ranging over several tens of meV. At each energy increment away from resonance, the de- 34 tuned photons can scatter inelastically, creating or annihilating phonons in the material, to obtain the resonance energy of the Mössbauer isotope (57Fe in our case). With some data reduction work [68, 69] the number of such inelastic, incoherent scattering events at each detuned energy yields the vibrational density of states. With NRIXS we obtain the partial density of states of the Mössbauer isotope only, since only these are excited when the detuned photons are put onto resonance through inelastic scattering. Comprehensive treatments of the NRIXS technique can be found in [52, 70]. 35 Chapter 3 Theoretical Foundations In this chapter we review common methods for solving the equations of quantum mechanics in materials. A basic understanding of quantum mechanics on the level of [71] is assumed. 3.1 Overview Materials are made of nuclei and electrons. The interaction between these basic entities determine all aspects of a materials behavior. The Hamiltonian for a system of interacting nuclei and electrons can be written [72], Ĥ = − X ~2 X 2 ZI 1 e2 X ZI ZJ e2 X ~2 X ∇2I ∇i + e2 + + − . 2me |ri − RI | 2 |ri − rj| 2 MI 2 |RI − RJ | i I,i i6=j I (3.1) I6=J Upper case indices refer to nuclear coordinates, while lower case indices are reserved for electrons. The terms in the Hamiltonian operator are defined as, 36 ~2 X 2 ∇i = electron kinetic energy 2me i e2 X I,i 2 X e 2 i6=j ZI = electron − nuclear Coulombic interaction |ri − RI | 1 = electron − electron Coulombic interaction |ri − rj| ~2 X ∇2I = kinetic energy of nuclei 2 MI I 2 X e 2 I6=J ZI ZJ = Coulomb interaction between nuclei. |RI − RJ | (3.2) If we could insert Eq. 3.1 into the expressions of quantum mechanics and solve the resulting equations, we could predict all properties of a material. This is a fulfillment of the reductionist dream of the ancient Greeks, who believed the entire phenomenal world to be reducible to the “atoms and the void” [73]. The fact that we are today the possessors of such a reductionistic theory of the phenomenal world is as awe inspiring as it is useless in practical terms. The equations of quantum mechanics are hopelessly complicated for a system of many interacting electrons and nuclei. Nobel laureate and the father of Density Functional Theory, Walter Kohn, has gone so far as to suggest that the many-electron wavefunction is “not a legitimate scientific concept” when the number of electrons exceeds about 103 [74]. The computational resources of the known Universe are quickly exhausted by the exponentially increasing resources required to store the many-body wavefunction as the number of interacting electrons grows. Fortunately, ingenious approximations have been devised that allow researchers today to cut the Gordian knot of quantum complexity and wring true predictive power out of the equations. In formulating these approximations, two general strategies have been pursued. In the first strategy, one retains much of the complexity in the wavefunction and concentrates on evaluating the equations very efficiently. An example of this strategy will be discussed later in this chapter when we briefly 37 review quantum Monte Carlo techniques. The second strategy attempts to reduce the intrinsic complexity of the quantum problem. A prime example of this strategy is density functional theory, to which we now turn. 3.2 Density Functional Theory Density functional theory reduces the complexity of the quantum many-body problem by decreasing the number of variables one must work with and then making clever approximations for the toughest terms in the simplified Hamiltonian. Our first task is to re-write the fundamental solid state Hamiltonian (Eq. 3.1) in a form more amenable to approximate solution. 3.2.1 Born-Oppenheimer Approximation The nuclei are much more massive than the electrons, and so move much more slowly. We expect atomic core (nuclear) velocities to not exceed 105 cm/s (the speed of sound), and the velocity of the electrons to be on the order of 108 cm/second [75]. This large disparity between nuclear and electronic motion means that the electrons experience the atomic cores as being essentially at rest at any given instant. The electronic system remains in its ground state throughout the atomic vibrations. So long as the quantum state of the system of electrons remains unchanged as the nuclei vibrate about their equilibrium positions, the nuclear motion can be decoupled from the electronic motion in the system. This decoupling is known as the Born-Oppenheimer, or adiabatic approximation. This approximation was first discussed with regards to molecular motion by Born and Oppenheimer in 1927 [76]. Reviews of the essential points of the Born-Oppenheimer approximation can be found in [77, 78, 79, 80, 81]. In practice, the Born-Oppenheimer approximation means that in studying the system of interacting electrons and nuclei we can freeze the nuclei in place. This gives, 38 ~2 X ∇2I =0 2 MI I 2 X e 2 I6=J ZI ZJ = constant. |RI − RJ | (3.3) The Hamiltonian of the electronic system now becomes, Ĥ = − X ~2 X 2 e2 X 1 ZI ∇i + + e2 + constant. 2me 2 |ri − rj| |ri − RI | i i6=j (3.4) I,i We can simplify our notation by writing this Hamiltonian as (following the notation of [72]), Ĥ = T̂ + V̂ext + V̂int + EII (3.5) where, T̂ = ~2 X 2 ∇i = electron kinetic energy 2me i V̂ext = e2 X I,i V̂int = EII = 2 X e 2 i6=j 2 X e 2 I6=J ZI = electron − nuclear Coulombic interaction |ri − RI | 1 = electron − electron Coulombic interaction |ri − rj| ZI ZJ = Coulomb interaction between nuclei. |RI − RJ | (3.6) The ext in V̂ext designates the potential the electrons feel from the nuclei as external, because the Born-Oppenheimer approximation removes the motion of the nuclear cores from the quantum problem. The electrons exist in the external potential the frozen nuclei create for them. This in contrast to the internal potential (Vint) the electrons feel from one another. 39 3.2.2 Hohenberg-Kohn Theorems: Charge Density as Fundamental Quantum Variable As discussed in the beginning of this chapter, the correlated wavefunction for a many-body system is a tremendously complex entity to work with. Ignoring spin and stipulating Ψ to be real, the many body wavefunction depends on 3N variables for a system of N interacting electrons. As N increases, the many-body wavefunction quickly becomes unmanageably complicated. However, quantum mechanics can be formulated in terms other than that of many-body wavefunctions. The important insight in this regard (for our purposes) began to crystallize in the fall of 1963 when two Americans, Walter Kohn and Pierre Hohenberg, came together at the École Normale Supérieure in Paris and began to look carefully at the fundamental importance of the electronic density in systems containing many electrons. An appreciation of the importance of the electron density distribution to the total energy of condensed matter systems developed as the two researchers considered the Hamiltonian of the CuZn alloy system [74]. Charge transfers occur between the Cu and Zn sites in disordered Cu-Zn alloys. These charge transfers lead to important electrostatic interaction energy contributions to the Hamiltonian of Cu-Zn alloys. This work led Hohenberg and Kohn to a careful consideration of the Thomas-Fermi theory, which attempts to express the electronic energy in the language of electronic density [82, 83]. The flaws in the Thomas-Fermi theory led Kohn to ask, “Is a complete, exact description of ground state electronic structure in terms of the electron density distribution possible in principle?” [74]. The answer came in the form of two theorems that would fundamentally reshape computational condensed matter physics. 40 First Hohenberg-Kohn Theorem: Vext(r) is (to within a constant) a unique functional of the electronic density distribution (ρ(r)) [84]. ρ =⇒ V̂ext =⇒ Ĥ =⇒ All P roperties The Hamiltonian (Ĥ) is determined by Vext(r), and so all of the properties of the many particle system follow from a specification of ρ(r). The second Hohenberg-Kohn theorem establishes a variational principle for the density functional. Second Hohenberg-Kohn Theorem: The lowest energy for the system is obtained if and only if the ground state density is used in the energy functional [84]. E [ρ0 ] ≤ E [ρ] Here we have indicated that the systems energy (E) is a functional1 of the electronic density distribution (ρ). Proofs of the Hohenberg-Kohn theorems are surprisingly simple. The first theorem is proved by postulating two (nontrivially) different external potentials, and showing that a contradiction can be derived (reductio ad absurdum) if both potentials result in the same electron density. The proof of the second theorem follows from the conventional variational principle of quantum mechanics (the one having to do with wavefunctions) [85]. With these simple proofs Hohenberg and Kohn shook the foundations of theoretical condensed matter physics and chemistry. The inconveniently complicated wavefunction no longer held center 1 A function maps one number to another number. A functional maps a function to a number [85]. The form of the function determines the value the functional gives. A functional is a ‘function of functions’. In this work (and generally) square brackets, G[f], are used for the argument of a functional (as opposed to the typical function notation f(x)). 41 stage. Turning to the electronic density as the fundamental quantum variable means moving from an entity (Ψ) that depends on 3N variables (N being the number of electrons in the system) to an entity that depends on only 3 variables (ρ). In conventional quantum mechanics we specify our system, which gives us the potential acting on the electrons (V ). This determines the Hamiltonian (H). In non-relativistic quantum mechanics the Schroödinger equation is used to find the wavefunction given the systems Hamiltonian. Once the wavefunction is determined the observables of the system are calculated by taking expectation values of the appropriate operators. The electronic density is one of these observables. Density functional theory turns this logical structure on its head. The electronic density, which is one observable among many in conventional quantum mechanics, is promoted to the foundational quantity of density functional theory. One begins with the electron density, which uniquely determines the potential acting on the electrons (this according to the Hohenberg-Kohn theorems). The potential then determines the Hamiltonian, which determines other observables (such as total energy E . . .). These logics are depicted in Fig. 3.1. Figure 3.1: Density functional theory promotes the electronic density, one among many observables in traditional quantum mechanics, to the primary quantum variable. Despite this great simplification in the foundations of materials theory, the very simplicity of the logical structure behind the Hohenberg-Kohn theorems is also the source of their greatest weakness. While asserting electronic density as fundamental, it fails to spell out how one might practically go about using this insight to do computational work in the real world. The theorems are not constructive. While the fundamental equation governing the evolution of the wavefunction is known, an equivalent equation for density is not known [86]. The next year in 1965, Kohn and Sham [87] would outline a practical method of putting the Hohenberg-Kohn Theorems to work. 42 3.2.3 Kohn-Sham Method 3.2.3.1 Kohn-Sham Hamiltonian The first step toward practical density functional theory calculations is to re-write the Hamiltonian in a form that facilitates approximate solution. We begin with the Hamiltonian after the BornOppenheimer approximation has been made, Ĥ = T̂ + V̂int + V̂ext + EII . (3.7) The general strategy in the Kohn-Sham method is to calculate as much of the Hamiltonian as possible exactly, and then make a clever and efficient approximation for those parts which cannot be easily calculated exactly. The kinetic energy of non-interacting electrons, the classical Coulomb electron-electron energy, and the classical Coulomb electron-nucleus energy can be calculated exactly. These terms are consequently included in our re-grouped Hamiltonian. We introduce them into our Hamiltonian by writing, Ĥ = T̂ + T̂0 − T̂0 + V̂int + V̂H − V̂H + V̂ext + EII (3.8) where T̂0 represents the kinetic energy of non-interacting electrons and V̂H is the Hartree potential for the classical Coulomb interaction between electrons and nuclei. The terms in the solid state Hamiltonian for which we do not have simple, closed-form expressions are grouped together into an exchange-correlation term,     V̂xc = T̂ − Tˆ0 + V̂int − V̂H . (3.9) Ĥ = T̂0 + V̂H + V̂ext + EII + V̂xc . (3.10) We can now write, 43 Equation 3.10 is the Kohn-Sham Hamiltonian. 3.2.3.2 Solving the Kohn-Sham Equations Iteratively As discussed earlier, the key difficulty in predicting the behavior of materials beginning from the full equations of quantum mechanics is the fiendishly complicated many-body wavefunction. In the Kohn-Sham method we skirt the many-body wavefunction and all of the difficulty it entails by making an audacious assumption. We assume that a system of non-interacting particles exists which gives the same ground state density as the fully interacting system of electrons present in our material. A priori, we have no reason to be confident that such a non-interacting system exists. Nevertheless, in practice, this assumption has proven successful beyond anyone’s reasonable expectation. The success of this approximation in practice may be gauged from the fact that the original paper by Kohn and Sham [87] laying out the method based on this assumption is the most cited Physical Review paper of all time [88]. It is even more amazing that the number of citations this paper receives has been growing every year for over forty years [88]. Assuming the existence of a non-interacting system of particles which gives the correct ground state density for the interacting system simplifies the many-body quantum problem to that of solving a series of independent-particle Schrödinger equations [89], Ĥ φi = i φi, 1≤i≤N (3.11) where Ĥ is the Kohn-Sham Hamiltonian, φi is the single-particle wavefunction for the i-th particle, N is the number of electrons in the system, and i are the eigenvalues of the Kohn-Sham Hamiltonian. The eigenvalue problem 3.11 is arrived at mathematically by using the variational principle to find the condition the single particle wavefunctions (orbitals) must satisfy for the energy functional to be at its minimum, given the orthonormality constraint [85]. After solving the independent-particle Schroödinger equation N times, the resulting orbitals are used to calculate the ground state electronic density of the interacting system, 44 ρ(r) = N X (φ∗i (r))(φ(r)). (3.12) i=1 To solve the Kohn-Sham equations the single-particle wavefunctions in Eq. 3.11 are expanded in a set of basis functions and the Kohn-Sham Hamiltonian is diagonalized. To get a feel for how this process works it will be helpful to look at the inner workings of an actual software package that solves the Kohn-Sham equations. To prepare ourselves for the chapters to come, we sketch the methods used in the popular Vienna ab initio simulation package (VASP) [90]. The Hohenberg-Kohn theorems assure us that the ground state electronic density uniquely determines the Hamiltonian, so the Kohn-Sham Hamiltonian is expressed in terms of ρ. This leads to a circular problem, in which the equations we must use to calculate ρ themselves depend on ρ. We address this conundrum by guessing an initial charge density and then iterating to self consistency. The VASP self-consistency cycle is depicted in Fig. 3.2. In the DFT calculations performed with the VASP package, the initial charge density guess is derived from a superposition of atomic charge densities. At the beginning of the calculation (step a in Fig. 3.2) the wavefunction arrays are populated with random numbers. The initial charge density is used to construct the Kohn-Sham Hamiltonian (step b). The eigenvalue problem is constructed by expanding the wavefunctions in a basis set of plane waves (step c). Applying Bloch’s theorem [75] we express the plane wave basis set as, 1 X CGnk ei(G+k)·r φn,k(r) √ V G (3.13) with k the wave vector, G a reciprocal lattice vector, n the band index, and V the unit cell volume. The band index, n, is used to differentiate the different independent eigenstates that arise for each k ([75]). The number of plane waves in a basis set is determined by specifying an energy cutoff (Ecut) in units of energy, so only those plane waves that satisfy, ~2 2 |G + k| < Ecut 2me (3.14) 45 Figure 3.2: Solving for the Kohn-Sham ground state through iteration to self-consistency. are included in the basis set. The Hamiltonian matrix is diagonalized at a number of k points in the Brillouin zone. The diagonalization of the Hamiltonian matrix (step d) gives the eigenvalues and eigenvectors according to the eigen decomposition theorem [91]. Matrix diagonalization is the major computational task that must be accomplished in real time during the course of a standard DFT calculation. Direct matrix diagonalization scales as the third power of the number of plane wave in the basis set. Realistic calculations often require hundreds to thousands of plane waves at each k-point. The number of k-points can easily number in the hundreds (which can often be reduced by taking symmetry into account). Given these significant computational requirements, direct matrix diagonalization quickly becomes intractable. In iterative matrix diagonalization, a simple process is repeatedly applied, pro- 46 ducing a series of solutions that converge to the desired solution. The iterative process ends when a pre-defined accuracy is achieved. Iterative matrix diagonalization methods are an order of magnitude faster than direct diagonalization and fit more naturally into the self-consistency cycle required by the circular density problem [90]. For these reasons iterative techniques are the standard in density functional theory calculations. To iteratively diagonalize the Hamiltonian matrix in the DFT problem, an expansion set is constructed that approximates the wavefunction. On each iteration, a residual vector is added to the expansion set to improve the approximation. The residual vector is a measure of the error between the wavefunction and the approximate representation [92]. Iterative diagonalization methods are differentiated by how the residual vector is added to the expansion set [92]. Common iterative diagonalization algorithms used in DFT calculations include the Blocked Davidson [93, 94, 95, 96], Conjugate Gradient [97, 98], and the Residual Minimization Method - Direct Inversion in the Iterative Subspace (RMM-DISS) [99, 100] algorithms. There are tradeoffs between these algorithms. For example, the RMM-DIIS algorithm converges on a solution the fastest, but can miss eigenvalues. The Blocked Davidson algorithm is somewhat slower, but is more stable than other choices. (The “Blocked” in “Blocked Davidson” implies that all orbitals are optimized simultaneously. Unblocked algorithms optimize bands one at time.) The algorithms can be combined for optimal performance [92]. After the wavefunctions are optimized in the iterative matrix diagonalization process, they are used to calculate a new charge density (step e). The old and the new charge density are then mixed (step f ), using, for example, the RMM-DIIS algorithm. This mixed charge density is used to construct a new Hamiltonian (step g) and calculate a new ground state energy. This new energy is compared to the input energy (step h). If the new energy is close enough to the input energy to satisfy the given consistency condition, the calculation ends and the total energy is returned. If the difference between the new and old energies is still too large, the mixed density is used as input for another round of iterative matrix diagonalization and density mixing. Now that we have some detailed understanding of the Kohn-Sham method, we can see that 47 employing the Kohn-Sham approximations essentially means that we are working with a mean field theory [86]. We approximate the hopelessly-complex fluctuating forces the electrons impose on one another with a steady, average force determined by the electronic density. This average force is used to solve for the behavior of each particle independently of the instantaneous behavior of the others. The fact that this simplification has often worked well enough in practice is fortunate indeed; but there is a price to be paid for this simplification. The seriousness of this price is determined by the degree of success we have in constructing an exchange-correlation term for our Kohn-Sham Hamiltonian. The true art in formulating a successful computational strategy based on density functional theory involves formulating a sufficiently accurate approximation for the exchangecorrelation term. 3.2.4 Exchange-Correlation Term The term exchange refers to the phenomenon of like-spin electrons (or Fermions in general) avoiding each other due to the antisymmetrization of the Fermionic wavefunction. Correlation refers to the coordination of movements among electrons in order to minimize Coulomb energy. In practice, the correlation portion of the Hamiltonian includes everything that is present in the exact solution but is unaccounted for in the exchange and classical Hamiltonian terms. The basic strategy of the Kohn-Sham formulation of density functional theory is to exactly express as much of the Hamiltonian as possible in terms of closed-form expressions. The remainder of the Hamiltonian, the most difficult parts, are then gathered into a single term, the exchangecorrelation term. The importance of this term cannot be overstated. If the exact form of the exchange-correlation term could be found, density functional theory would be an exact ground state many-body theory. As exciting as this may sound, an exact closed form expression for the exchangecorrelation functional likely does not exist. Nevertheless, highly successful approximations to the exchange-correlation functional have been discovered, and are the crucial ingredient that has allowed density functional theory to become a workhorse for the computational materials community. 48 3.2.4.1 Local Density Approximation The exchange-correlation energy of a homogeneous electron gas was calculated successfully by Ceperley and Alder in 1980 [101] using the quantum Monte Carlo technique (of which more will be said later in this chapter). Viewed on the scale of a unit cell, a typical material bears little resemblance to a homogeneous electron gas (see Fig. 3.3), making Ceperley and Adler’s success seemingly irrelevant in real materials. The idea behind the local density approximation is to divide the unit cell into very small volume elements and use Ceperley and Alder’s result within each small volume element. Within sufficiently small volume regions, the electron density will be approximately uniform and Ceperley and Alder’s results (and their interpolation to other volumes) may be relevant. Adding up the contribution to the exchange and correlation energy (xc ) from each volume element d3r then gives the exchange-correlation term for the entire material (Exc ), Exc = Z ρ(r) homog (ρ(r)) d3 r, xc (3.15) where homog is the exchange-correlation energy of a homogeneous electron-gas at density ρ. The xc exchange-correlation potential is given as the functional derivative of the exchange-correlation energy (Exc ) with respect to the local electronic density, V̂xc = δExc [ρ] . δρ (3.16) The local density approximation is generalized to magnetic systems by setting, Exc = Z ρ(r) homog (ρ ↑ (r), ρ ↓ (r)) d3r, xc (3.17) where ρ ↑ (r) and ρ ↓ (r) are the densities of spin-up and spin-down electrons, respectively. In this local spin density approximation (LSDA), the different spin densities are treated as two independent variables. The correlation term accounts for interactions between electrons of opposite spin [102]. 49 Figure 3.3: Electronic isosurface (mesh) around atomic cores (solid spheres) in a unit cell shows complex electron behavior in a nonhomogeneous system. If the variation in electronic density is sufficiently slow, the exchange and correlation from a nonhomogeneous system can be captured by local quantities. This idea forms the basis of the local density approximation. 3.2.4.2 Generalized Gradient Approximation The source of the exchange-correlation energy in the local (spin) density approximation is both the source of the L(S)DA’s greatest strength and its greatest weakness. The fact that the exchangecorrelation energy is derived from very accurate calculations on a uniform electron gas means that both the exchange and the correlation terms satisfy physically important normalization conditions [102]. This helps to explain the broad success of the L(S)DA despite the simple nature of the underlying approximation. Conversely, being based on calculations on a uniform electron gas, the LDA cannot be completely accurate for non-homogeneous systems. In non-homogeneous systems the exchange-correlation term depends both on the electron density at a point r and also on its variation in the immediate vicinity of r. The LDA neglects these gradients and so is most accurate for systems with slowly varying electronic density. The generalized gradient approximation (GGA) attempts to improve on the LDA by including gradient terms to account for the variation of the density around r [72], 50 Exc = Z ρ(r) xc (ρ(r), ∇ρ) d3r. (3.18) The GGA formulation is not uniquely defined (in contrast to the LDA), and several formulations exist. The GGA is generalized to magnetically-polarized systems in a manner analogous to the LDA. The GGA succeeds in some areas where the LDA has been found lacking; but it is not a universal improvement over the local approximation. In general, the GGA improves on the underestimation of the exchange energy and the overestimation of the correlation energy endemic in the LDA. The opposite nature of these corrections results in inconsistent improvements as one treats different materials, however [102]. Careful testing and comparison to measured quantities for a given material are needed before a reasoned choice between the LDA and GGA can be made. 3.2.5 The Projector Augmented Wave Method In Section 3.2.3.2 the expansion of the wavefunctions in a plane wave basis set was discussed. Expanding electronic wavefunctions in plane wave basis sets is an efficient strategy for valence electrons due to their free electron-like nature. Nearer to the nucleus, the Coulomb potential from the positively charged ions becomes significantly stronger, and the wavefunctions of the core electrons are consequently much more oscillatory in nature. It is not efficient to expand the core state wavefunctions in plane wave basis sets. This unfortunate situation is helped by the fact that the core electrons are largely unaffected by the bonding behavior of the atoms, which primarily involves the valence electrons. This means that the behavior of condensed matter systems can largely be derived from the behavior of the valence electrons alone. The core electrons do not play a significant role apart from generating, together with nuclei, the potential in which the valence electrons move. This means that if we could replace the core electrons with an alternative system that produces identical effects outside the core, we could greatly simplify the electronic structure problem. This is precisely what is attempted by the use of pseudopotentials in electronic structure calculations. The scattering of a particle can be formulated in terms of the phase shift in the scattered 51 particle’s wavefunction caused by the scattering potential. The key theoretical point on which all pseudopotential methods rest is the fact that the scattered wavefunction outside the core region is indistinguishable for phase shifts modulo 2nπ (n an integer) [72]. Consequently, we are free to replace the significant nuclear Coulomb potential by a much weaker potential, so long as it produces the same phase shift modulo 2nπ on scattering. The formulation of pseudopotentials is not unique [72]. This being the case, the goal in formulating a pseudopotential is to find one that is useful in the chosen computational scheme and that accurately reproduces the scattering phase shift modulo 2nπ over the target energy range [72]. While the formulation of pseudopotentials giving accurate results in condensed matter systems has been the foundation for many of the advances in computational materials science over the last several years, the elimination of core electron states does not come without some consequence. It is sometimes the case that core states do matter. In these cases, the retention of the core electronic states is essential. The price we pay for keeping all of the electrons in the calculation is a significant increase in the computational resources we must employ to treat a given system accurately. A nice compromise has been formulated recently that allows for the efficiency of pseudopotential methods while retaining the core electron wavefunctions. This breakthrough is the Projector Augmented Wave (PAW) method [103, 104]. The PAW method begins with the expansion of the wavefunctions in a plane wave basis set, as with traditional pseudopotential-plane wave methods. The new idea in the PAW method is to use a set of projectors on radial grids at the atom centers to re-introduce the core wavefunctions into the calculation in what essentially amounts to a frozencore approximation. This procedure retains access to the full all-electron wavefunction and charge density while maintaining the efficiency of the pseudopotential methods. The key innovation in the PAW method is the use of radial grids together with an angular momentum representation to evaluate the atomic core terms [104]. The regular grids used in the atomic core regions in earlier pseudopotential methods made the retention of the all-electron wavefunction too computationally costly to implement. The increased computational efficiency stemming from the use of radial grids in the PAW method allows one to work directly with the all-electron wavefunctions 52 and potentials [104]. To understand how the PAW method works, let us first introduce a smooth wavefunction, |Ψ̃i (employing Dirac notation [71] for convenience). We expand the smooth (pseudo) wavefunction in spherical harmonics (partial waves) within muffin-tin spheres around each atomic site, E X E ψ̃ = cm ψ̃m (3.19) m with the corresponding all-electron wavefunction, |ψi = X cm |ψm i . (3.20) m E Now we let |ψi and ψ̃ be related by the linear transformation τ , |ψi = τ ψ̃ E (3.21) We can now write (substituting 3.19 into 3.20), E X n Eo |ψi = ψ̃ + cm |ψm i − ψ̃m (3.22) m τ is linear, so the coefficients in the expansion must be given by a projection carried out in each sphere, for a projector set p̃ E D cm = p̃m ψ̃m . (3.23) Substituting 3.23 into 3.22, and regrouping, we obtain, E Xn Eo D E |ψi = ψ̃ + |ψm i − ψ̃m p̃m ψ̃ . m Comparing 3.21 to 3.24 we see immediately that, (3.24) 53 τ = 1+ Eo Xn |ψm i − ψ̃m hp̃m | . (3.25) m Equation 3.24 shows us that the linear transformation decomposes the wavefunction into three terms within the augmentation spheres. An exact onsite atomic core wavefunction is added to the smooth (pseudo) wavefunction. The smooth atomic site wavefunction is then subtracted from the overall wavefunction (see Fig. 3.4). Figure 3.4: The projector augmented wave method (see text for explanation). (wf denotes wavefunction.) In the PAW method the linear transformation between the all-electron and pseudo (smooth) wavefunction is also applied to the charge densities and the total energy functional. When this is done, it is found that an analogous decomposition into three terms also holds for each term in the energy functional, as well as for the charge densities [104]. 3.3 First-Principles Lattice Dynamics The ability to simulate lattice vibrations from first-principles allows one to predict thermodynamic properties central to materials behavior. This includes the ability to extend 0 K DFT results to finite temperatures. An example of this type of calculation is given in Appendix B. There exist several approaches to first-principles lattice dynamics. These include the linear response method [105, 106, 107], the frozen-phonon method [108] and the small-displacement method [109, 110]. For 54 the calculation of full dispersion curves and the vibrational density of states, the linear response and small-displacement methods are the most efficient. In this work we will concentrate on the small-displacement method. While both the frozen phonon and the small-displacement method also go by the name direct method, we will always refer to the small-displacement method when using the term direct method. A brief review of lattice dynamics will facilitate an understanding of the direct method. Extended reviews of lattice dynamics can be found in [81, 78]. Atoms in a crystal vibrate about their equilibrium lattice positions, even at 0 K. We let V(R(n,µ), . . . R(m,υ)) denote the ground state many-body potential of a system of atoms at lattice sites R(n, µ), where n indexes the primitive unit cell, and µ indexes the atomic site within the unit cell. We can then express the instantaneous position of an atom as it vibrates about its equilibrium position R as, r = R + u, (3.26) with u denoting the momentary displacement of the atom from its equilibrium position. Expanding V about u we have, V (R(n, µ), ...R(m, υ)...) = V0 + X (u(n, µ) − 0) n,µ, + ∂V ∂u(n, µ) 0 X 1 ∂2V (u(n, µ) − 0) (u(m, υ) − 0) . 2 n,µ,m,υ ∂u(n, µ)∂u(m, υ) 0 (3.27) The many-body potential, V, exists as a definite quantity only if the motion of the electrons and the motion of the atomic cores can be decoupled. We will assume that the adiabatic approximation (see section 3.2.1) is valid throughout this derivation. The second term is identical to zero as there are no forces on the atoms at equilibrium, ∂V = 0. ∂u(n, µ) 0 (3.28) 55 We define the elements of the force constant matrix as, Φi,j (n, µ, m, υ) = ∂2V . ∂u(n, µ)∂u(m, υ) 0 (3.29) Φi,j (n, µ, m, υ) is the negative of the force in the i direction on the atom at (n, µ) when the atom at (m, υ) is displaced in the j direction. There will be one 3x3 force constant matrix for each pair of atoms that interacts in the material. Working in the harmonic approximation, we express the potential as, V (R(n, µ), ...R(m, υ)...) = V0 + 1 X u(n, µ) u(m, υ) Φi,j (n, µ, m, υ). 2 n,µ,m,υ (3.30) The equations of motion for the lattice are obtained from Newton’s Second Law, M ü(n, µ) = −∂V , ∂u(n, µ) (3.31) and so we have, M ü(n, µ) = − X Φ(n, µ, m, υ)u(m, υ). (3.32) m,υ, The equations of motion are solved with a plane wave solution, 1 u(n, µ) = p (k, j) ei(k·R(n,µ)−ω(k,j)t). Mµ (3.33) Here Mµ is the atomic mass, (k, j) is the wave amplitude, k is the wave vector, j indexes the mode, ω(k, j) is the phonon frequency for wave vector k and mode j. R(n, µ) is the real space position of atom (n, µ). Substituting 3.33 into 3.32 and rearranging, we define the dynamical matrix as, X 1 D(k; µ, υ) = p Φ(0, µ; m, µ)exp(−2πik · [R(0, µ) − R(m, υ)]), Mµ Mυ m (3.34) 56 summing over all atoms in the crystal. We can now express the equations as an eigenvalue problem using matrix notation, ω2 (k, j)(k, j) = D(k)(k, j). (3.35) Here D(k) is a 3r x 3r matrix, with r the number of atoms in the basis. (k,j) is a column vector of size 3r. The diagonalization of the dynamical matrix at each k yields 3r eigenvectors () and 3r eigenvalues (ω). The small-displacement method assumes that Hooke’s law holds for atoms in a crystal; namely that forces on displaced atoms scale linearly with displacement amplitude. Starting from Eq. 3.30 this assumption gives the expression, Fi(n, µ) = − X ∂V Φi,j (n, µ, m, υ)uj (m, υ), =− ∂u(n, µ) (3.36) m,υ,j relating the displacement (u(m, υ)) of an atom (m, υ) from its equilibrium position to the forces (F(n, µ)) generated on all other atoms due to this displacement. In the small-displacement method, the determination of the force constant matrix (Φi,j (n, µ, m, υ)) is the goal, as it will allow us to calculate the lattice dynamical behavior. The Hellmann-Feynman theorem [111, 112] allows us to extract forces (Fi (n, µ)) from the displacement of atoms (uj (m, υ)) in a supercell through a simple fixed-point energy calculation on supercells with single atoms displaced from equilibrium. The theorem states that the force, F, on an atom at position R is given by, 57 ∂E ∂R E ∂ D =− Ψ Ĥ Ψ ∂R * + ∂ Ĥ =− Ψ Ψ . ∂R F=− (3.37) This by the extremal property of the exact wavefunction [72]. The Hellmann-Feynman theorem allows us to efficiently extract the forces on all other atoms due to the displacement of one atom in a supercell. We just need to calculate the ground state electronic density in the displaced cell. Having carried out the calculation of the ground state electronic density and extracting the Hellmann-Feynman forces (Fi(n, µ)) for a given atomic displacement (uj (m, υ)), the elements of the force constant matrices (Φi,j (n, µ, m, υ)) can be extracted. In this thesis we employ the lattice dynamics software package PHONON [113] in conjunction with the density functional theory package VASP [90] to obtain lattice dynamical behavior from first-principles. PHONON finds the optimal (in the least-squares sense) force constant elements using the singular value decomposition method [113, 114]. Each symmetry inequivalent atom is displaced in each symmetry inequivalent direction. A supercell is created for each atomic displacement and the Hellmann-Feynman forces are computed for each atom in the distorted supercell. To avoid error due to the periodic nature of the boundary conditions, the supercell must be large enough so that the force constants fall off by at least three orders of magnitude at the edge of the supercell compared to the center of the supercell centered on an atom [113]. After their extraction through the fitting process to the Hellmann-Feynman forces, the force constant matrix is used to construct the dynamical matrix (Eq. 3.34). The diagonalization of the dynamical matrix at a series of wave vectors (k) along straight lines in reciprocal space gives the phonon dispersion curves for the crystal. The diagonalization of the dynamical matrix at a set of randomly distributed (selected by Monte Carlo sampling) k vectors distributed homogeneously over 58 the first Brillouin zone yields the phonon density of states. The relation between polarization vectors ((k)) and the direction of wave propagation (k) divide phonon modes into those of longitudinal and transverse nature. In an isotropic medium, for a given k, it is always possible to choose solutions such that one of the modes is polarized along the direction of wave propagation (kk), while two are polarized perpendicular to the direction of wave propagation (⊥k). We call the branch with kk the longitudinal branch. The branches with ⊥k are the transverse branches. An additional classification of phonon modes derives from the behavior of a phonon branch as k vanishes at the origin of the Brillouin Zone. Acoustic branches are those that show a vanishing frequency as k goes to the origin. There are always exactly 3 acoustic branches to match the 3 dimensions of space we experience. 3r-3 of the branches are optical branches, which maintain finite frequencies of vibration as k goes to 0. These optic branches tend to be flatter due to the fact that in these modes the intracellular interactions are stronger than the interactions between cells. 3.4 Quantum Monte Carlo In the introduction to this chapter we alluded to two broad strategies used to solve the quantum solid-state problem. Density functional theory exemplified the strategy of reducing the number of quantum variables through approximation. Quantum Monte Carlo methods are an example of an alternative strategy. They attempt to retain the full quantum equations and instead concentrate on efficiency of evaluation. Working with the full many-body problem means performing numerical integrations involving many-dimensional wavefunctions. Random sampling is far and away the most efficient method of carrying these evaluations out [86]. The use of random sampling is the string that ties all Monte Carlo techniques together. Despite an every increasing efficiency, quantum Monte Carlo techniques remain more computationally expensive than DFT techniques. If DFT can provide the accuracy needed to solve a problem, then it is generally the most efficient choice of ab initio techniques. Quantum Monte Carlo techniques are reserved for those systems on which the approximations of density functional theory fail catastrophically. 59 In this section we will sketch two quantum Monte Carlo techniques that are of particular, and growing, importance. These are the variational and diffusion quantum Monte Carlo techniques. We will then, through an example calculation, show the power of diffusion Monte Carlo in treating systems involving weak interactions. 3.4.1 Variational Quantum Monte Carlo In the variational quantum Monte Carlo (VMC) technique we simply guess the form of the manybody wavefunction and then evaluate the expectation value of the Hamiltonian operator using Monte Carlo techniques to obtain the energy of the system [86]. The guessed wavefunction generally depends on a set of parameters that can be varied to minimize the energy and optimize the wavefunctions form. In variational quantum Monte Carlo calculations one typically uses the Slater-Jastrow wavefunction [115], Ψ(R) = D eJ . (3.38) D is a Slater determinant of electronic orbitals, D= φ1(x1 ) ··· φ1 (xN ) .. . .. .. . φN (x1) . (3.39) · · · φN (xN ) where N is the number electrons in the system. The Jastrow factor, eJ , explicitly includes terms that account for electron correlation, and is optimized in the course of a calculation. The more parameters we are able to optimize in our trial wavefunction, the more accurate our final result is likely to be. This can lead to complicated expressions for the trial wavefunction. In general, we do not expect to capture all of the many-body effects in variational quantum Monte Carlo calculations; 80-90% of a system’s correlation energy is the maximum we can generally expect to capture [86]. This is generally not good enough to justify the computational cost involved in 60 VMC calculations. VMC is most useful as a means of generating an optimized wavefunction with which to begin the diffusion quantum Monte Carlo technique. 3.4.2 Diffusion Quantum Monte Carlo If we make the substitution τ = it in the Schrödinger equation we obtain the imaginary-time Schrödinger equation,   1 ∂Ψ(R, τ ) . − ∇2 + V (R)) Ψ(R, τ ) = − 2 ∂τ (3.40) Now if we ignore the potential term we obtain the diffusion equation, 1 ∂Ψ(R, τ ) − ∇2Ψ(R, τ ) = − . 2 ∂τ (3.41) Ignoring the kinetic term, we obtain the rate equation, V (R))Ψ(R, τ ) = − ∂Ψ(R, τ ) . ∂τ (3.42) Seen in this way, the imaginary time Schrödinger equation describes walkers in 3N dimensional space diffusing and dying or reproducing at a pace set by the value of the potential V (R) [116]. The wavefunction (Ψ) is represented by an ensemble of walkers [86]. In time the walker density comes to be proportional to the ground state wavefunction for the system, and accurate expectation values can then be calculated. This is the essence of the diffusion Monte Carlo (DMC) algorithm [117, 118]. In theory, the DMC method can capture all of the correlation effects in a many-body system. The major limitation to achieving this goal in a fermionic system is that, left to itself, the DMC algorithm always converges to the ground state bosonic wavefunction [86]. For a system of Fermions (such as electrons) the fixed-node approximation must be imposed. For Fermions, the wavefunction must be an antisymmetric function of the Fermion coordinates. This means that the wavefunction must change signs. In making the fixed-node approximation [119] we impose a nodal surface on the wavefunction. Across the nodes of the nodal surface the wavefunction changes signs. Once the nodal 61 surface is constructed, the DMC algorithm is carried out within each region bounded by nodes. The guessed structure of the nodal surface is a fundamental barrier to complete accuracy in the DMC method; improving on the fixed-node approximation is an active area of research in the quantum Monte Carlo community. In practice in a DMC calculation, one begins with a trial wavefunction, often the Slater-Jastrow wavefunction as in VMC. A common method is to populate the Slater determinant with singleparticle orbitals generated in a preliminary density functional theory calculation on the system. The nodal surface is also taken from this initial DFT calculation. Variational quantum Monte Carlo calculations are carried out to optimize the trial wavefunction and to create an initial configuration of walkers for the DMC method to work with. The DMC algorithm is then carried out on the population of walkers until the desired accuracy is achieved. 3.4.3 Quantum Monte Carlo Example: C2 H4 The study of systems involving van der Waals interactions are difficult to carry out computationally due to the origin of van der Waal forces in the correlated motions of many electrons and the weak nature of the interaction. Traditional DFT methods treat the weak van der Waals interactions notoriously poorly. The unsuitableness of traditional DFT methods for van der Waals interactions can be traced to the local approximations these DFT techniques are based on. A true many-body technique, such as DMC, is needed to accurately capture the subtle electronic effects underlying van der Waal systems. The van der Waal forces between hydrogen atoms and carbon materials is a topic of current interest in the hydrogen storage community. Small molecular systems can give insight into behaviors likely to play important roles in larger carbon systems. To this end, we have computed the van der Waals interactions between an H atom and ethylene (C2H4 ) using ground state quantum Monte Carlo methods as implemented in the CASINO software package [120, 121, 86]. Fixed-node [119] diffusion quantum Monte Carlo [117, 118] calculations using a trial wavefunction of Slater-Jastrow form were performed. The orbitals in the Slater wavefunction were generated from DFT-PBE cal- 62 Figure 3.5: Potential energy of placing a hydrogen atom at various distances above the surface of a C2 H4 molecule, as shown in Fig. 3.6. The dots are the results of DFT calculations. The crosses are results from the diffusion quantum Monte Carlo calculations, which handles better the long-range Coulomb correlations of electrons in C2H4 and H. This more accurate QMC calculation predicts a much stronger potential for sorption at 3.6 Å. culations performed with the ABINIT package [122]. The Jastrow correlation function allows for explicit inclusion of long range correlation effects in the trial wavefunction. The Jastrow factor included electron-electron, electron-nucleus, and electron-electron-nucleus terms. The inclusion of electron-electron terms in the trial wavefunction is vital in order to capture the long-ranged correlations responsible for dispersion (van der Waals) forces. The free parameters in the Jastrow factor were optimized in preliminary variational quantum Monte Carlo calculations. The optimized many-electron wavefunction from VMC was then used in the diffusion quantum Monte Carlo calculations. The imaginary-time Schrödinger equation was used to evolve an ensemble of electronic configurations, improving the ground state estimate beyond that which is possible with DFT. 63 3 Figure 3.6: Electron density contour (0.23 q/Å ) of a hydrogen atom above an ethylene molecule at a height of 2.25 Å (left) and 2.0 Å (right). These are density functional theory calculations, and correspond to the situation just above and below the barrier at 2.1 Å shown in Fig. 3.5. The results in Fig. 3.5 show that the inclusion of long-range electron-electron correlations cause a much stronger van der Waals attraction of hydrogen on ethylene at a distance of 3.6 Å than do standard DFT calculations (also shown in Fig. 3.5). The barrier at 2.1 Å prevents strong bonding to the Carbon atoms. The fact that the barrier is larger in the QMC result than it is in the DFT calculation has potentially important implications for the storage of hydrogen on carbon structures for use in mobile applications, for example. Figure 3.6 shows that once this barrier is overcome, electron transfer between the hydrogen and the ethylene becomes important as a chemical bond is formed. For the vertical parabolic potential shown in the inset of Fig. 3.5, we calculate a vibrational frequency for a H atom of 34.15 meV. The validity of this quantum Monte Carlo result could be probed experimentally using inelastic neutron scattering. 3.5 Conclusion With the steady beat of Moore’s Law behind it, computational materials physics will continue from its already prominent position to assume ever greater importance in the toolbox of the materials scientist in the years ahead. The greatest utility of these techniques are found when they are used in close conjunction with experimental work; guiding the experiments forward and verify and rationalizing the experimental results. The brief review of techniques in this chapter gives a taste of the 64 methods and approximations that are being used today to tackle some of the toughest computational problems in science. The field is a vibrant one, and many innovative techniques are on the horizon that will make tomorrow a bright one indeed for the computational researcher. 65 Chapter 4 Pressure-Induced Invar Behavior in Pd3Fe The Invar effect has remained a topic of interest since Guillaume discovered the low thermal expansion of Fe-Ni alloys in 1897 [4]. Iron-rich fcc-based Invar alloys are characterized by the existence of nearly degenerate states with differing magnetic moments and volumes [123] (and refs. therein). Decreasing the Fe concentration stabilizes the magnetically ordered state, and suppresses thermal Invar behavior. Recently it was discovered [1] that some Fe-Ni alloys with normal thermal expansion properties at ambient pressure exhibit Invar behavior at high pressure. Alloys of Fe-Pd with compositions around 30% Pd have long been known for Invar behavior at ambient pressure [2]. In addition to Invar behavior, Fe-Pd exhibits several other anomalies in mechanical and magnetic behavior such as a martensitic transformation in Fe-rich samples [124], giant magnetic moment formation [125], noncollinear magnetic states [126], and anisotropy in the spin-wave dispersions in Pd-rich alloys [127]. In this chapter we report pressure-induced Invar behavior in Pd-rich Pd3 Fe with the ordered L12 structure, shown in Fig. 4.1. At ambient pressure, Pd3 Fe is ferromagnetically ordered with a Curie temperature (Tc ) measured between 500 K [128] and 540 K [129], and magnetic moments at room temperature of 2.73±0.13 µB and 0.51±0.05 µB at the Fe and Pd sites, respectively [130]. We report that the ferromagnetic ground state is destabilized with pressure, and a magnetic collapse begins around 10 GPa, accompanied by an increase in bulk modulus at pressures after the magnetic collapse. We also find that at a pressure near the start of the magnetic collapse, ordered Pd3 Fe 66 exhibits a near-zero thermal expansion from 300 to 523 K. Figure 4.1: L12 crystal structure for Pd3 Fe. 4.1 Experimental Methods and Results 57 Samples of Pd57 Fe of 95.38% isotopic 3 Fe were prepared by arc-melting Pd of 99.95% purity and enrichment. There was no measurable mass loss or visible surface oxidation after melting the 100 mg ingots. The material was cold rolled to a thickness of 25 µm. The L12 chemical long-range order was obtained by annealing at 700◦C for 16 hours, followed by annealing at 600◦C for 24 hours, and a subsequent cooling to 20◦C over 2 hours. X-ray diffraction confirmed the L12 structure with a rolling texture consistent with (011)[211̄]. Room temperature x-ray diffraction (XRD) measurements were performed at pressures up to 33 GPa at beamline X17C of the National Synchrotron Light Source (NSLS) at Brookhaven National Laboratory (BNL). Energy-dispersive x-ray diffraction (EDXD) data were acquired with a polychromatic x-ray beam and a Ge detector at a fixed Bragg angle of 12◦ . A Merrill-Bassett diamond-anvil cell [48] with 500 µm culets and a 250 µm chamber was used with silicone oil as the pressure medium, with pressure calibration by the fluorescence of ruby technique [131]. No crystallographic phase transition was observed under pressure. 67 Figure 4.2: Volume-pressure data obtained from synchrotron x-ray diffraction measurements (symbols) and the fit of the data to the Weiss-like equation of state [20] (line). Figure 4.2 shows volume vs. pressure data for ordered Pd3 Fe from the EDXD measurements. A significant volume collapse was found between approximately 10 and 15 GPa, and a reduced slope in the volume-pressure data was found above 15 GPa. The EDXD data were fit successfully to the “Weiss-like equation of state” [20], V (P ) =  1 1 + e(P −PT )/∆P   + 1−    0 V0,HS 1 + B0,HS 1 1 + e(P −PT )/∆P   P B0,HS  0 V0,LS 1 + B0,LS 0 −(1/B0,HS ) P B0,LS   0 −(1/B0,LS )  , (4.1) 0 with volume V corresponding to pressure P, B0 the zero-pressure isothermal bulk modulus, B0 the pressure derivative of the bulk modulus, V0 the ground state volume, PT the pressure of the transition between the high-spin (HS) and low-spin (LS) regions, and ∆P the pressure range over which the volume collapse takes place. 68 The high-spin (HS) and low-spin (LS) regions are modeled with two separate Murnaghan [51] equations of state,  −(1/B00 ) 0 P V (P ) = V0 1 + B0 . B0 (4.2) The region around the volume collapse is modeled as a weighed sum of the high-pressure and lowpressure equations of state, with a thermodynamic occupancy factor n(P), n(P ) = 1 1 + e(P −PT )/∆P , (4.3) based on the assumption that the energies of the electronic energy levels have a linear dependence on pressure. Further details on the Weiss-like equation of state are given in Section 2.2.1. The pressure-volume data divides into three regimes, labeled in Fig. 4.2. In the low-pressure region (region I in Fig. 4.2) up to 10 GPa the fit to the equation of state gave a zero-pressure bulk modulus (B 0 ) of 229±2 GPa and a ground state lattice parameter (a0 ) of 3.849±0.002Å, the latter result in good agreement with previous studies [132]. The high-pressure region from 16.5 to 33 GPa (region III) gave a B 0 of 286±3 GPa and an a0 of 3.792±0.005Å. The pressure derivative of the bulk modulus was fixed to the typical metallic value of 4.0 in all cases. In the volume-collapse region (region II) from 10 GPa to 16.5 GPa, fitting to the Weiss-like equation of state gave a transition pressure between the low-pressure and high-pressure regions of 12.1±0.1 GPa and a pressure range for the transition of 0.9±0.1 GPa. Nuclear forward scattering (NFS) spectrometry offers a direct measure of the magnetic state under pressure, and was performed at 300 K at beamline 16ID-D at the Advanced Photon Source (APS) at the Argonne National Laboratory (ANL) using the same diamond anvil cell [48]. Figure 5.1 shows NFS spectra from ordered Pd57 3 Fe at pressures up to 25.3 GPa. The transition region of Fig. 5.1 follows well the transition region shown in Fig. 4.2. Quantum beats, expected from a magnetically-ordered material, are prominent in the NFS spectra at lower pressures. The quantum beats decrease significantly in amplitude at 12.3 GPa, and vanish at the highest pressures. They 69 Figure 4.3: Nuclear forward scattering spectra from ordered Pd57 3 Fe showing an abrupt decrease of the hyperfine magnetic fields in the material around 12 GPa. do not change significantly in their time modulation, however. This is consistent with a first-order phase transition, where pressure replaces a high magnetization state with a low magnetization state. In another EDXD experiment, the Pd57 3 Fe sample was heated under pressure in a diamond anvil cell with a resistive heating furnace (D’Anvils Ltd.) (Figure 4.4) at beamline X17C of the NSLS. Au powder was used as an internal pressure calibrator inside the cell during the heating. Measurements were taken at a pressure of 7 GPa (±0.1 GPa at the lower temperatures), heating to a maximum temperature of 673 K. Figure 4.5 shows the results of these measurements. The Au expands at all temperatures while the Pd3 Fe shows practically no thermal expansion up to 523 K. (Above 600 K there was a more rapid loss of pressure in the cell. We expect this would reverse some of the Invar 70 transition leading to an artificially large thermal expansion at the highest temperatures of Fig. 4.5.) Figure 4.4: DAC furnace on NSLS X17C beamline (left) with cutaway of furnace and DAC (right). Cutaway graphic created by Lisa Mauger. 4.2 Density Functional Theory Calculations We performed floating moment density functional theory (DFT) calculations in the local spin-density approximation (LSDA) [101] for the electronic exchange and correlation potential using the Vienna Ab-initio Simulation Package (VASP) [90]. The interaction between the electrons and ions was treated using the projector augmented-wave (PAW) method [133], treating the semicore 2p states in Fe and 3p states in Pd as valence. For the plane wave basis set, an energy cutoff of 800 eV was used, the large value facilitating accurate calculation of the stress tensor. The k-space summation was performed using a 17 × 17 × 17 Monkhorst–Pack grid [134], giving 165 k-points in the irreducible portion of the Brillouin zone. The partial occupancies were determined according to the first–order Methfessel–Paxton [135] method with a 0.1 eV smearing width, giving a fictitious electronic entropy term not exceeding 1 meV/atom. Total energies were converged to better than 1 meV/atom. Calculations were performed for a ferromagnetic state (FM), a state of ferrimagnetic order (FR), in which the spin on the central Fe atom in a 32 atom supercell is reversed with respect to the other 71 Figure 4.5: Volume-temperature data for Pd57 3 Fe (diamonds) and Au (squares) obtained from externally heating the diamond anvil cell at an approximately constant pressure of 7 GPa. The Pd57 3 Fe shows practically no thermal expansion up to 523 K. The solid and dashed lines are guides for the eye. spins, a low-spin state (LS), and three antiferromagnetic states, which we denote AFM-I, AFM-II, and AFM-III. The Fe site spin directions associated with the magnetically ordered states are shown in Fig. 4.6. In the LS state studied the moments at the Fe and Pd sites couple ferromagetically with their combined magnetic moment in the unit cell never exceeding 0.01 µB . In the AFM-I state, Fe moments order in alternating ferromagnetic sheets in the (110) plane. The AFM-II state has Fe moments ordering ferromagnetically in alternating (100) planes, while in the AFM-III state they order ferromagnetically in alternating (111) planes. Different energy cutoffs for the plane wave basis sets, and different k-space grids were used for the various magnetic structures owing to the different supercells. In all cases, total energies were converged with respect to both k-points and energy cutoff to better than 1 meV/atom. For each magnetic state, energy versus volume data were obtained by first relaxing the geometry according to the conjugate-gradient algorithm, then performing a series of single-point energy 72 Figure 4.6: Spin directions at Fe sites in the calculated magnetically ordered states in Pd3 Fe. calculations about the equilibrium volume. The resulting energy-volume data were fitted to the Murnaghan equation of state [51], giving parameters listed in Table B.1. Table 4.1: Parameters obtained from DFT calculations of energy versus volume, fitted to the Murnaghan equation of state. FM AFM-I AFM-II AFM-III FR LS Fit Range(Å3 /atom) 12.33 ≤ V ≤ 15.07 12.21 ≤ V ≤ 14.92 12.48 ≤ V ≤ 14.66 12.54 ≤ V ≤ 15.00 12.27 ≤ V ≤ 15.00 11.84 ≤ V ≤ 14.48 B0 (GPa) 216.81 216.39 217.79 218.19 218.06 242.59 B00 5.12 5.03 5.02 5.15 5.06 5.11 a0 (Å) 3.802 3.798 3.800 3.800 3.801 3.760 This procedure gave a FM ground state with a lattice parameter of 3.802 Å at 0 GPa and 0 K with magnetic moments of 3.16 and 0.33 µB within spheres of radius 1.486 Å around the Fe and Pd sites (similar to earlier spin-polarized DFT results [136, 13]). The calculated lattice parameter is slightly smaller than the experimental result, as is expected for the LSDA, while the magnetic moments are in good agreement with saturation moments [137]. Figure 4.7 shows that over a large volume range from 9 to 15 Å3 /atom, the FM, AFM-I, AFM-II, 73 and FR states are all extremely close in energy, never separated by more than 9 meV/atom, and have similar elastic properties. The AFM-III structure is within 37 meV/atom in this same volume range. With increasing pressure the FM state is destabilized, leading to a transition to the AFM-I state at a=3.690 Å. Charge is transferred to the Fe site as the lattice parameter decreases, as shown in Fig. 4.8. In Fig. 4.9 we see that the magnetic moments in the AFM-I state decrease gradually with decreasing volume, disappearing completely by a=3.148 Å. If only the FM and LS states are taken into account, a transition to the LS state takes place at a=3.355 Å, where the FM moment on the Fe site drops discontinuously from 1.39 µB to 0.008 µB . Such behavior would be qualitatively similar to that predicted for disordered Fe70 Pt30 [138], with a larger stability range for the AFM state in Pd3 Fe. Figure 4.7: Energy-volume data from DFT LSDA calculations on Pd3 Fe. Inset shows a close-up of the low-pressure region. 4.3 Discussion While stable only at very high pressures, the elastic properties of the calculated LS state show similarities to the high-pressure state in the experimental studies. Figure 4.10 shows an excellent match 74 Figure 4.8: Charge at the Fe site versus the lattice constant of Pd3 Fe from DFT calculations. Antiferromagnetic refers to the AFM-I structure (see text). The dashed line is a guide to the eye. between the pressure-volume data from EDXD measurements and the DFT calculated equations of state for the ferromagnetic and low-spin states. There is a significantly higher B 0 and lower V 0 in the LS state. These changes upon transition to the LS state are similar to those of the transition found in the EDXD data and the NFS data of Figures 4.2 and 5.1, but the calculated transitions occur at very different pressures. The EDXD measurements show a volume collapse beginning around 10.1 GPa, or V /V 0=0.96, ending by 14.8 GPa, or V /V 0 =0.91, similar to that of the FM to AFM-I transition found in our DFT calculations (V/V0 =0.91). The DFT calculations do not give a stable LS state until a much higher compression (V/V0 =0.57). In the LSDA some discrepancy between theoretical and experimental transition pressures is expected, and a higher pressure is typical, but this discrepancy is larger than expected. We discuss two possible interpretations. The simplest is that the LSDA is greatly underestimating the relative stability of the LS magnetic state, and a high-spin (HS) to LS transition occurs around 75 Figure 4.9: Magnetic moments at the Fe site versus the lattice constant of Pd3 Fe from DFT calculations. Antiferromagnetic refers to the AFM-I structure (see text). The solid line is a guide to the eye. V /V 0 =0.96 instead of the calculated V /V 0=0.57. This view is supported by the similarities in the changes in B 0 and V 0 between the HS and LS states in the calculations and experiments. Further support is given by the excellent fit to the EDXD data with the Weiss-like equation of state, which is based on the assumption of the Weiss 2-γ model [6] of a transition between a HS state at low pressures to a LS state at high pressures. An alternative explanation is possible. There will also be a collapse of the 57Fe hyperfine magnetic field if the moments fluctuate at a rate of 109 Hz or higher. At ambient pressure the magnetic ordering temperature in Pd3 Fe is ∼500 K. Pressure is known to have a strong effect on the value of Tc in Invar alloys [5, 139] (and refs. therein). Increasing the pressure above 10 GPa may reduce Tc to below room temperature in Pd3 Fe, where a spin disordering transition can occur. Nevertheless, the Fe moments in Pd3 Fe are well localized [13], and the Fe atoms likely retain most of their moments in the 76 Figure 4.10: EDXD pressure-volume data (diamonds) together with DFT-calculated equation of state data for the ferromagnetic and low-spin states. paramagnetic state. This paramagnetic state is to be distinguished from a truly non-magnetic state in which no local moments exist in the material [140]. The magnitude of the discrepancy between the experimental and theoretical transition pressures leads us to favor this alternative explanation. Nevertheless, a significant change in the local magnetic moment magnitude could accompany this transition, as discussed for Fe-Pt and Fe-Pd [141]. A quantitative account of the spontaneous volume magnetostriction (or the fractional change in equilibrium volumes through the magnetic transition, ωs), ωs = V (F erromagnetic) − V (P aramagnetic) V (P aramagnetic) (4.4) in Invar has been given for Fe-Pt and Fe-Pd [141] using the disordered local moment method (DLM) [34]. These studies identified the necessary condition for Invar behavior as having an alloy at the edge of the transition from strong to weak ferromagnetism, i.e., the Fermi level (E F ) must be at 77 the top of the majority 3d band in the ferromagnetic state. Thermal excitations cause a relatively large decrease in the local Fe moments, and large magnetovolume effects around the transition temperature. Our calculated results in Fig. 4.11 show that alloying with Pd causes the top of the majority 3d band to move to energies below E F , stabilizing the strong ferromagnetic state and suppressing the Invar effect, in agreement with [141]. On the other hand, our calculations in Fig. 4.12 show that applying external pressure gives the opposite trend, as the antibonding states at the top of the majority 3d band move up in energy towards E F . Pressure therefore counteracts the band-filling effect of Pd, bringing an anomalously large spontaneous volume magnetostriction to Pd3 Fe. By tuning the position of the top of the 3d band with respect to E F , pressure-induced Invar behavior resembles classical Invar behavior that is controlled by chemical composition. Figure 4.11: Majority-spin electronic density of states at the Fe site for alloys from the Fe-Pd system at ambient pressure. 78 Figure 4.12: Pd3 Fe majority-spin electronic density of states at the Fe site at various volume compressions (Fe3 Pd is shown for reference). 4.4 Conclusion Our theoretical work in this chapter shows that the DLM theory of classical Invar behavior extends successfully to the pressure-induced case. Electronic structure calculations point to a unification of the classical Invar and pressure-induced Invar pictures. Our results on Pd3 Fe also answer, for pressure-induced Invar behavior, a question that has been of historical importance for classical Invar alloys. Invar behavior was originally discovered in chemically disordered Fe64Ni36. Subsequently, a long debate ensued as to the role chemical disorder does or does not play in Invar behavior. Our results on Pd3 Fe show for the first time that chemical disorder is not a pre-requisite for pressureinduced Invar behavior. Pd3 Fe is the first instance of pressure-induced Invar behavior in a chemically ordered material. The work related in this chapter brings pressure-induced Invar behavior outside the Fe-Ni family for the first time. There is no reason to believe it will be the last. Pressure-induced Invar behavior 79 is likely a common phenomenon in alloy systems containing a member of the Invar family. Pressure, then, becomes an additional means by which researchers can tune to an Invar anomaly. By searching for additional cases of pressure-induced Invar behavior and understanding its microscopic underpinnings in each case, researchers will bring a powerful new light to bear on the century old question of the origin of the Invar effect. 80 Chapter 5 Iron atom dynamics across the pressure-induced Invar transition in Pd3Fe 5.1 Introduction The Fe-Pd alloy system exhibits interesting properties owing to a subtle interplay between magnetism and mechanical behavior. These properties, which are highly sensitive to chemical composition, include the shape-memory effect [142], martensitic instability [143], and vanishing thermal expansion (Invar effect) [144]. Interatomic forces are central to these phenomena, and studies of lattice dynamics can reveal some of the underlying reasons for the relationships between mechanical behavior and magnetism. Stirling, et al. [145], measured phonon dispersions in ordered Pd3 Fe using inelastic neutron scattering. Fitting the results to a Born-von Kármán model, they found weaker nearest-neighbor Pd-Fe force constants, compared to the nearest-neighbor Pd-Pd force constants. This difference was attributed to the size mismatch between Fe and Pd atoms and the distance dependence of the d-electron interaction. Yue, et al. [146], measured the 57Fe phonon partial density of states (PDOS) in L12 -ordered Pd3 Fe and Pt3 Fe using nuclear resonant inelastic x-ray scattering (NRIXS), and reported effects from the mass difference of Pd and Pt atoms. The optical modes in Pd3 Fe involving Fe atoms were broadened in energy, compared to those in Pt3 Fe. Ghosh [147] performed a theoretical study of the lattice dynamics in ordered Pd3 Fe using density functional perturbation 81 theory (DFPT) [148]. The calculated phonon DOS and force constants were in reasonable agreement with the earlier experimental studies. The control of pressure is a powerful means to alter both lattice dynamics and the magnetic state of materials. Pressure has long been used in studies of Invar materials [149, 150, 151, 152], but recently a new phenomenon was discovered at high pressures in Fe-Ni alloys away from the classical Invar composition [1]. The application of pressure transformed the non-Invar alloys Fe0.55Ni0.45 and Fe0.20Ni0.80 into Invar alloys with zero thermal expansion. Several of the authors recently found this pressure-induced Invar phenomenon in Pd3 Fe, which is far-removed from the narrow composition range of the Invar alloy near Fe0.70Pd0.30. Under a pressure of 7 GPa the alloy Pd3 Fe exhibits negligible thermal expansion from 300 K to approximately 523 K [153]. The same study showed that the volume collapse at pressures from 10 and 15 GPa was accompanied by a magnetic transition from the ferromagnetic (FM) ground state to a low-spin (LS) state. Density functional theory (DFT) calculations showed how pressure moved the majority-spin antibonding t2g levels across the Fermi surface, accounting for the magnetic transition. The calculations did not give a reasonable prediction for the pressure of the magnetic transition, however, and this was attributed to an inadequate accounting for the spin structure in the high-pressure LS phase of Pd3 Fe. Here we report results from an investigation on the lattice dynamics across the pressure-induced Invar transition in Pd3 Fe. The NRIXS method was used to measure the 57Fe PDOS for samples in diamond anvil cells at pressures up to 35 GPa. The DFT calculations of the PDOS were successful for the FM state at low pressure, but were less successful for the low-spin state at high pressures. For the NRIXS spectra measured at high pressures, the highest energies of the phonon spectra showed features consistent with the DFT predictions for FM material, suggesting there may be short-range magnetic order in the high-pressure phase of Pd3 Fe. 82 5.2 Methods 5.2.1 Nuclear Forward Scattering Nuclear forward scattering measurements were performed at beamline 16ID-D at the Advanced Photon Source at the Argonne National Laboratory. The preparation of Pd3 57Fe by the arc-melting of small ingots was described previously [153]. A Merrill-Bassett, Tel-Aviv-type diamond anvil cell (DAC) [48] was used to generate pressure on the sample. The sample was compressed in a sample chamber created by a Re gasket (∼100 µm diameter hole) and two opposing diamond culets of ∼350 µm diameter. A 4:1 methanol-ethanol mixture was used as the pressure medium in the sample chamber. For low temperatures, the DAC was cooled with a He-flow cryostat (Cryo Industries of America, Inc.). Two thermocouples measured temperatures throughout the low-temperature runs. An in-line Raman system was used for in situ pressure calibration at cryogenic temperatures by means of ruby fluorescence [49], and this method was also used for pressure determination in cells at ambient temperature. Figure 5.1 shows quantum beats in the NFS time spectra, indicative of a magnetically-ordered material at lower pressures and ambient temperature. As the pressure is increased to 16.6 GPa at 300 K, the quantum beats largely disappear, indicating the collapse of magnetic order under pressure. These new measurements on Pd3Fe are consistent with our earlier results [153]. At a fixed pressure of 16.6 GPa, the quantum beats reappear upon reducing the temperature, although their periodicity is shorter, owing to a larger hyperfine magnetic field at cryogenic temperatures. As the temperature is increased at 16.6 GPa, the quantum beats gradually fade, disappearing by 100 K. The absence of quantum beats is also consistent with the presence of magnetic moments if they fluctuate faster than approximately 1 ns. This time scale is long compared to the vibrational periods of the phonons measured by nuclear resonant x-ray scattering. We also know that temperature suppresses the quantum beats at a temperature of 500 K as Pd3 Fe passes through its Curie transition [128], so we can give a rough picture of the P − T phase boundary for the Invar transition in Pd3 Fe. Around ambient temperature and below, the 83 temperature-dependence is weak, with the transition occurring at approximately 16.0 GPa at 300 K, and approximately 16.6 GPa at 90 K. The temperature dependence must be stronger at higher temperatures, however, reaching 0 GPa at 500 K. These numbers should be considered as approximate, however. In one sample there was hysteresis in the magnetic transition, with the LS state persisting to low pressures after it formed at high pressure. Figure 5.1: 57Fe nuclear forward scattering spectra from L12 -ordered Pd3 57Fe. Quantum beats, indicative of the low-pressure magnetically-ordered phase with moments stable longer than 1 ns, collapse above 15.5 GPa. Quantum beats reappear on lowering the temperature to 5.1 K, and persist up to 75 K at 16.6 GPa. 84 5.2.2 Nuclear Resonant Inelastic X-Ray Scattering Using similar procedures as for the NFS samples, samples of Pd3 57 Fe were prepared by arc-melting Pd and 57 Fe of 95% isotopic enrichment. The sample was cold rolled into a thin foil of approximately 30 µm thickness, and annealed in vacuum to induce L12 long-range order. X-ray diffraction measurements confirmed prominent superlattice diffraction peaks characteristic of the L12 structure. Nuclear resonant inelastic x-ray scattering (NRIXS) spectra were measured at the 3-ID undulator beamline of the Advanced Photon Source at the Argonne National Laboratory. Final monochromization to a 1 meV bandwidth was performed with a high-resolution silicon (4 0 0)(10 6 4) monochromater [64]. Samples were loaded into a piston-cylinder type diamond anvil cell [154] for measurements to 35 GPa. The cell was sealed with Be gaskets, 100 µm in thickness. Silicone oil was used as the pressure medium. Pressure was measured using the ruby fluorescence technique [49] before and after the acquisition of spectra at each pressure. Spectra were collected by tuning the x-ray energy stepwise in the range of ±80 meV around the 57Fe Mössbauer nuclear resonance energy. Twelve scans were performed at each pressure, each scan approximately one hour, and the successful scans were summed. Three avalanche photodiodes were used to count the delayed photons from nuclear de-excitation. The inelastic scattering intensity from Pd3 57Fe was converted into a 57Fe phonon partial density of states (PDOS) with the PHOENIX software package [68]. Results are presented in Fig. 5.2. With increasing pressure, the PDOS spectrum splits. The high-energy peak shifts with pressure rapidly to higher energies. 5.2.3 Born-von Kármán Fitting The experimental 57Fe phonon PDOS curves were fit to a Born-von Kármán model [155, 156] with first- and second-nearest neighbor (1NN, 2NN) interactions. A set of interatomic force constants was taken from interpolated DFT results, and used as the initial guess for the iterative procedure to obtain force constants. The calculated PDOS was convolved with the experimental instrument resolution function when comparing to the NRIXS spectra. The force constants were optimized to converge on the measured NRIXS PDOS result, using a differential evolution algorithm [157] 85 Figure 5.2: 57Fe phonon PDOS from NRIXS measurements on L12 -ordered Pd3 Fe. The Fe PDOS from a fit of the 0 GPa results to a Born-von Kármán (BvK) model is also shown. implemented by the Mystic software package [158]. The Pd phonon PDOS is not measurable with NRIXS, so the Pd-Pd force constants are poorly constrained when fitting the Born-von Kármán model to 57Fe phonon PDOS curves. The 1NN and 2NN Pd-Pd force constants were therefore fixed to the Pd-Pd force constants from the DFT calculations, interpolated to the appropriate experimental volumes. Table B.1 shows that these DFT Pd-Pd force constants match reasonably well those from inelastic neutron scattering measurements [145] at ambient pressure. The resulting force constant matrices were decomposed into longitudinal and transverse components in two steps [159]. First, the force constant matrices were diagonalized, yielding three eigenvalues. Next, the force constant matrices were projected onto the bond directions. The eigenvalue matching the projection is the longitudinal force constant. The other two eigenvalues were averaged to give an averaged transverse force constant. This approach is reasonable so long as the forces are primarily axial [159], as is typical of fcc metals. 86 Table 5.1: Force constants projected along bond directions (longitudinal) and average perpendicular (transverse) at V/V0=1.0. [145] 1NN Pd-Pd (expt.) 1NN Pd-Pd (calc.) [145] 2NN Pd-Pd (expt.) 2NN Pd-Pd (calc.) 5.2.4 Longitudinal (N/m) 57.28 59.52 5.21 3.46 Transverse (N/m) -5.448 -5.00 1.945 -0.01 Density Functional Theory Calculations The Vienna Ab-initio Simulation Package (VASP) [90] was used for calculations on L12-ordered Pd3 Fe with density functional theory in the local spin-density approximation (LSDA) [101] for the electronic exchange and correlation potential, using the projector augmented-wave (PAW) method [133]. In the calculations the semicore 3p states of Pd and 2p states of Fe were treated as valence. A unit cell was constructed and relaxed, followed by a series of single-point energy calculations around the equilibrium volume. The calculated results were fit to a Murnaghan equation of state [51], and the ground state volume, bulk modulus and the pressure derivative of the bulk modulus were obtained, as described previously [153]. To study the vibrational properties of L12 -ordered Pd3 Fe, a series of 32-atom supercells were constructed. In the supercells, each symmetry-inequivalent atom was displaced from its equilibrium position in each symmetry-inequivalent direction. The Hellmann-Feynman forces induced on all other atoms were determined from a total energy calculation. An energy cutoff of 330 eV was used for the planewave basis set, together with a 6×6×6 Monkhorst-Pack grid. The first-order MethfesselPaxton method with a smearing width of 0.1 eV was used for k-space summations. Energies were converged to better than 2 meV/atom with respect to both planewave cutoff and k-point grid . The resultant Hellmann-Feynman forces were used to generate the dynamical matrix using the direct method [110] as implemented in the PHONON software package [113]. The diagonalization of the dynamical matrix at a series of points in reciprocal space gave the phonon dispersion curves and the phonon density of states (DOS). Supercells for ferromagnetic (FM), antiferromagnetic (AFM), low-spin (LS), and nonmagnetic (NM) states were constructed. In the AFM state, the Fe moments 87 were ordered in alternating ferromagnetic sheets in the (110) plane. Moments at the Fe and Pd sites were coupled ferromagetically in the LS state with the combined magnetic moment in the unit cell never exceeding 0.1 µB. For calculations on materials under pressure, volume was constrained in the supercell calculations to better reveal the effect of the magnetic state on vibrational properties. To account for instrument resolution broadening and to allow for direct comparison with experimental results, the calculated vibrational spectra were broadened by a convolution with a Gaussian function with a standard deviation of 1 meV. 5.3 Results and Discussion 5.3.1 Average Phonon Energy The first moment, hEi, of the Fe PDOS is shown in Fig. 5.3 as a function of decreasing volume (with respect to the volume at a pressure of 0 GPa, V0 ). The DFT results fall on two curves, a high-moment curve with stiffer phonon energies for the ferromagnetic and antiferromagnetic states, and a low-moment curve with softer phonons for the LS and NM states. In Fig. 5.3a the highmoment DFT solutions reproduce the NRIXS results well up to V/V0 =0.95 (12 GPa). Beyond V/V0 =0.95 the low-moment DFT solutions are more successful. The NRIXS hEi curve makes a transition between the high-moment and low-moment DFT curves, and this occurs at the pressures of the magnetic transition from the high-moment low-moment states found in the nuclear forward scattering measurements shown in Fig. 5.1. The slow stiffening of the 57Fe PDOS from V/V0 =0.95 (12 GPa) to V/V0 =0.89 (24 GPa) accompanies the magnetic transition under pressure in Pd3 Fe. In Fig. 5.3 we see that the FM and AFM states have similar average vibrational properties. This is also true for the LS and NM calculations. 5.3.2 Force Constants The Fe-Pd and Fe-Fe force constants for first- and second-nearest neighbors (1NN and 2NN) obtained by fitting the NRIXS curves to a Born-von Kármán model are presented in Fig. 5.4. Our 1NN 88 Figure 5.3: Measured and calculated average vibrational frequencies as functions of reduced volume for Fe (upper) and Pd (lower) phonon PDOS. The dashed lines are trend lines for the ferromagnetic and low-spin data, drawn as guides to the eye. Fe-Pd longitudinal force constants at ambient pressure (20.9 N/m) are in reasonable agreement with Stirling’s results (23.6 N/m) [145]. The large 1NN Fe-Pd longitudinal force constants stiffen significantly with pressure, but show a change in slope between 12 and 24 GPa. This pattern follows the trend of hEi of the NRIXS Fe PDOS shown in Fig. 5.3a. The smaller 0 GPa 1NN Fe-Pd transverse force constants (0.74 N/m) have opposite sign and are somewhat larger than Stirling’s (-0.4 N/m) [145]. The Fe-Fe 2NN longitudinal force constants are much smaller than the Fe-Pd 1NN longitudinal force constants. The Fe-Fe 2NN transverse force constants are large and negative. The longitudinal force constants calculated by DFT methods for the FM and LS states were 89 Figure 5.4: Fe force constants as a function of reduced volume calculated from fits of the NRIXS results to a Born-von Kármán model. in reasonable agreement with the results of Fig. 5.4 for Born-von Kármán fitting to the NRIXS spectra, with the 2NN longitudinal force constant being much smaller. Unfortunately, the methods provide different transverse force constants, although these are smaller. Nevertheless, we expect the longitudinal force constants to dominate the lattice dynamics. 5.3.3 Electronic Structure Figure 5.5 shows DFT results for the electronic DOS at the Fe site for L12 Pd3Fe, decomposed into states of t2g and eg character. (The Pd site is less interesting for the magnetic transition, since it has a small moment in the ferromagnetic state.) The ferromagnetic state was stable at larger 90 volumes, such as V/V0=1.0 (i.e., the relaxed structure at 0 GPa) and V/V0 =0.86. At V/V0 =0.85 the calculation has been performed in the low-moment state with vanishing magnetic moments at the Fe and Pd sites. These features of the electronic structure are similar to those proposed by Kaspar, et al. [38] and Entel, et al. [8] for classical Invar alloys. For fcc transition metal alloys, the t2g orbitals form strong dd bonds between 1NN atoms, owing to their large charge density in the [110] direction. This causes the t2g electronic DOS to split into a high-energy sub-band with anti-bonding character and a lowenergy sub-band with primarily bonding character. Both are occupied by majority spin electrons in the ferromagnetic state. (States with eg character primarily have charge density between the more distant 2NN atoms along the [100] direction. They are less affected by pressure.) With decreasing lattice parameter, there is an increased splitting of the bonding and anti-bonding t2g states. The occupied majority-spin anti-bonding t2g states move closer to the Fermi energy, as shown in Fig. 5.5, and the pressure-induced Invar transition occurs as they become energetically unfavorable and their electrons are transferred into minority spin states. 5.3.4 Lattice Dynamics Across the Invar Transition Figure 5.6a shows ferromagnetic Fe PDOS curves measured by NRIXS and calculated by DFT for V/V0 =0.99. The NRIXS result is broader, but agreement between experiment and theory is good for this magnetic state. Figure 5.6b compares the measured Fe PDOS to the calculated FM and LS curves at V/V0 =0.89. The NFS data of Fig. 5.1 show that Pd3 Fe is in a nonmagnetic state at a compression of V/V0 =0.89, so it might seem best to compare the experimental curve to the LS DFT curve. Although the low-energy part of the measured Fe PDOS matches the LS result reasonably well, in the LS curve the high-energy peak is too weak in intensity and energy. Figure 5.6b shows that turning on ferromagnetism in the DFT calculation causes a stiffening and significant increase in the relative intensity of the calculated high-energy peak in the Fe PDOS. The average of the FM and LS PDOS calculations seems to match better the Fe PDOS at high energies. The electronic states of t2g character involve dd electron bonding along 1NN directions. The 91 Figure 5.5: Electronic density of states at the Fe site decomposed into t2g and eg symmetry at various cell compressions. EF = 0 eV. 1NN forces are dominant, so the changes in t2g states across the Invar transition may dominate the changes in the lattice dynamics. With the collapse of the FM order comes a significant increase in the number of electronic states at the Fermi level, N (EF ), including more states with t2g character. More states at the Fermi level allows more efficient screening of the atom displacements in phonons, softening them. This is seen as a major change in the Fe PDOS for the calculated FM to LS transition (Fig. 5.6), but this change is too large compared to experiment. Again, we suggest that short-range magnetic order in the LS state may allow for stiffer 1NN force constants between Fe and Pd. Figure 5.7 shows phonon dispersion curves along high-symmetry directions, generated by the 92 Figure 5.6: Fe PDOS from NRIXS measurements and first-principles calculations. V0 for the NRIXS curves refers to the ground state volume determined by XRD in [153]. V0 for the DFT curves refers to the ferromagnetic ground state volume as determined by LSDA calculations [153]. force constants obtained from Born-von Kármán model fits to the NRIXS spectra (Fig. 5.7a), and by DFT methods for the LS and FM states (Figs. 5.7b,c). To the right of the dispersions are shown Fe PDOS curves, obtained as the sum of the mean-squared vibrational amplitudes of Fe atoms in all modes. The magnitudes of the red dots in Fig. 5.7 are proportional to the square of the vibrational amplitudes of the Fe atoms in selected modes. Figure 5.7 shows that the most important contribution to the high-frequency peak in the Fe PDOS originates from the high-energy modes around the Brillouin zone boundary (X and M points). The phonon modes between 30 and 35 meV around the X-point in Fig. 5.7a contribute substantially to the high-frequency peak in the Fe PDOS. Figure 5.8 depicts the Fe atomic displacements 93 Figure 5.7: Dispersion curves and Fe phonon PDOS for the (a) NRIXS measured, (b) LS and (c) FM DFT calculated states. See text for explanation. associated with the mode at 33.2 meV at the X-point in Fig. 5.7a, and this mode is marked with horizontal blue lines in Fig. 5.7. This mode involves planes of Fe atoms moving directly against one another within the unit cell. This mode and its neighbors dominate the high-frequency peak in the Fe PDOS in the NRIXS results at V/V0=0.89 (Fig. 5.7a). The DFT LS curve (Fig. 5.7b) predicts it to be too small, whereas calculations for the ferromagnetic state show a peak at 34 meV that is somewhat stronger than measured (Fig. 5.7c). This difference can be traced to a stiffening of this mode in the FM state, but also to a 70% decrease in the squared vibrational amplitude of the Fe atom motions in the LS state. This comparison of DFT and NRIXS results suggest that short-range magnetic order persists in 94 Figure 5.8: Two snapshots of the atomic displacements in a Pd3 Fe unit cell associated with the high-frequency mode at the X-point of the Brillouin zone (see text). Arrows indicate the motion of the Fe atoms. Pd3 Fe beyond the pressure-induced collapse of long-range ferromagnetic order. The phonon modes involving localized Fe atom movements, such as shown in Fig. 5.8, may be sensitive to short-range ferromagnetic order that exists in the LS state. The contribution of these modes to the Fe PDOS lies between the LS and FM DFT predictions. On the other hand, the long wavelength modes may be better modeled by the LS DFT predictions. There is other evidence for short-range magnetic order in Pd3 Fe. Non-collinear magnetic states with correlation lengths between three [160] and five [161, 162] times the lattice parameter have been reported in disordered Pd-Fe alloys with Fe concentrations around 15 at.%. Similar results, with correlation lengths extending over the length of a few lattice parameters, have been reported for Pt-Fe [161], Pt-Co, Pd-Co [163], Pd-Cr [164], and Pd-Mn [165] alloys. 5.4 Conclusions The 57Fe phonon partial density of states (PDOS) in L12 ordered Pd3 Fe was studied at high pressures using nuclear resonant inelastic x-ray scattering (NRIXS) measurements and density functional theory (DFT) calculations. On the average, the Fe PDOS stiffens with pressure more slowly between 12 and 24 GPa, owing to the pressure-induced Invar transition in Pd3 Fe. First-nearest neighbor 95 (1NN) Fe-Pd longitudinal force constants, obtained by fitting a Born-von Karman model to the experimental 57Fe PDOS, showed a relative softening between 12 and 24 GPa. The same general behavior was accounted for by the DFT calculations on ferromagnetic (FM) and low-spin (LS) states. At low pressures, the Fe PDOS calculated by DFT methods for the FM state matched well the measured 57Fe PDOS. At high pressures, for the higher-energy part of the Fe PDOS, the shape calculated for the LS state differs significantly from the measured spectrum, which seems to lie between the calculated shapes for the LS and ferromagnetic states. This discrepancy for higher-energy phonons is consistent with some short-range magnetic order in the LS state at high pressures. Nevertheless, the LS calculations give better agreement with the average phonon energy. Both the fits to experimental spectra and the DFT calculations show that Fe-Pd 1NN longitudinal force constants are stiffened by ferromagnetism. We interpret the change in lattice dynamics across the magnetic transition with the same changes in the electronic states of t2g character near the Fermi level that are important for the pressureinduced Invar transition itself. The energies of the antibonding t2g states increase relatively rapidly with pressure, and the occupied t2g states of majority spin become increasingly unfavorable compared to other electronic states. These states become depopulated at pressures of the Invar transition, and there is an increase in the number of t2g states at the Fermi level. This change at the Fermi level allows for better screening of the charge disturbances of atom vibrations, and a softening of the interatomic force constants. Nevertheless, although the DFT calculations work well for predicting the Fe PDOS in the ferromagnetic state, the predicted change in the Fe PDOS with the Invar transition is too large, especially for phonons of short wavevectors. This may be evidence of shortrange magnetic order in the LS state. 96 Chapter 6 Concluding Remarks The complexities and richness of the Invar state have fascinated and perplexed researchers for over a century. Our understanding of Invar behavior deepens and new questions arise with every advance in experimental and theoretical technique. As we move into the heart of this new century much work remains ahead on the Invar problem. In Chapter 4 we presented our discovery of pressure-induced Invar behavior in Pd3 Fe. This discovery marks the first instance of a pressure-induced Invar anomaly outside of the Fe-Ni family, and the first time a chemically ordered material has been found to display pressure-induced Invar behavior. Pd3 Fe exhibits a large volume collapse at room temperature beginning at 10 GPa, and collapse of hyperfine magnetic fields at this same pressure as revealed by nuclear forward scattering measurements. The high-pressure phase above 15 GPa is characterized by an increased bulk modulus in both the x-ray diffraction measurements and the density functional theory calculations. Under a pressure of 7 GPa, Pd3 Fe displays vanishing thermal expansion between room temperature and 523 K. Throughout the volume range studied, first-principles calculations reveal states of differing magnetic order competing closely in energy. The calculated ferromagnetic ground state is destabilized at high pressures, giving way first to an antiferromagnetic state, followed by a low-spin state at very high pressures as local moments decrease continuously with increasing pressure. When only the ferromagnetic and low-spin states are accounted for, the local magnetic moments vanish discontinuously at high pressure. Accompanying the decrease in magnetic moment magnitude at high 97 pressures is a transfer of charge from the larger Pd atoms to the smaller Fe atoms in Pd3 Fe. Interpreted in terms of the main Invar theories presented in Chapter 1, the experimental and computational results presented in Chapter 4 give a mixed picture. The interpretation of the pressureinduced Invar anomaly in Pd3 Fe in terms of the Weiss 2-γ model is supported by several of the results. First, the high-pressure x-ray diffraction data were fit successfully to the Weiss-like equation of state, which is based on the Weiss 2-γ model. Second, the nuclear forward scattering results point to a first-order magnetic phase transition with pressure, further supporting the Weiss 2-γ picture. The experimental results show a similar, if somewhat more extreme, stiffening in the low-spin state to the density functional theory results. This similarity gives a third piece of support for the Weiss 2-γ model. The measured and calculated magnetic transition pressures are very different, though, and this fact casts serious doubt on any simple interpretation of the experimental results in terms of traditional Weiss 2-γ ideas. In particular, the density functional theory results make doubtful the proposition that the simple density functional theory representation of the low-spin state as a state containing no local moments is adequate in the case of Pd3 Fe at high pressures. A viable alternative to the Weiss 2-γ picture is found in the disordered local moment model, which incorporates spin-disordering into its representation of the low-spin state. Our density functional theory calculations in the ferromagnetic state show that pressure counteracts the band-filling effect of Pd. Beginning from the classical Invar composition, Fe3 Pd, increasing Pd content at ambient pressure stabilizes the strong ferromagnetic state, moving the top of the majority 3d band away from the Fermi level. Our density functional theory calculations show that with increasing pressure the majority-spin t2g (see Fig. 6.1) anti-bonding electronic states move closer to the Fermi level. This restores the Pd3 Fe electronic structure to the condition responsible for Invar behavior in classical Invar alloy Fe3 Pd at ambient pressure. With this insight, the first step is taken toward a unification of the theoretical descriptions underlying classical and pressure-induced Invar behavior. The inelastic x-ray measurements and ab initio calculations presented in Chapter 5 extend our experimental and theoretical insights concerning pressure-induced Invar behavior in Pd3 Fe to the 98 Figure 6.1: (upper) t2g and eg orbitals in Pd3 Fe. (lower) The contribution of the t2g and eg electronic states to the majority-spin electronic density of states at the Fe site in Pd3 Fe. Ef =0 eV. dynamics of the crystal lattice. Once again, the results did not fit neatly into any one theoretical Invar framework. The measured average energy of the Fe phonon partial density of states stiffens with pressure, with the stiffening showing a slowing across the magnetic transition pressure. This trend is mirrored in the first nearest-neighbor Fe-Pd longitudinal force constants obtained from fitting the measured nuclear resonant inelastic x-ray scattering spectra to a Born-von Kámán model. The ferromagnetic density functional theory phonon results match the experimental results for average energy well at low pressures, while the low-spin density functional theory phonon results match the experimental results well at high pressures. These results again support an interpretation in terms of a simple 99 Weiss 2-γ picture. An examination of the electronic density of states reveals that the same movement of the t2g majority-spin antibonding states towards the Fermi level with increasing pressure that is responsible for pressure-induced Invar behavior explains the change in average phonon energy and interatomic forces that accompany the pressure-induced Invar effect. Through the transition to the low-spin state, t2g majority spin states move across the Fermi level, increasing the number of electronic states at the Fermi level, resulting in more effective screening of the atomic motions associated with lattice vibrations. This insight resonates with both the Weiss 2-γ model and the disordered local moment picture. The situation becomes more complicated when a detailed comparison of the shapes of the calculated and measured high-pressure Fe partial phonon density of states is made. The shape of the low-spin density functional theory result around the high-energy peak in the Fe phonon partial density of states disagrees with the measured result at high pressure. The situation is partially remedied by imposing ferromagnetic order at high pressures in the calculations, suggesting the possibility that short range magnetic correlations persist well beyond the magnetic transition pressure in Pd3 Fe. The insights into the pressure-induced Invar state yielded through modern synchrotron methods and density functional theory techniques are substantial. Nevertheless, an honest evaluation of the results presented in this thesis suggests that the power of the computational tools employed were not sufficient to model the high-pressure magnetic state in Pd3 Fe with complete accuracy. Evidence of their failure can be found in their inability to predict an accurate magnetic transition pressure in Chapter 4, and their failure to capture the behavior of the high-energy optical phonon modes at high pressures in Chapter 5. Both of these failures can be traced to an oversimplified representation of the high-pressure, low-spin magnetic state in Pd3 Fe. The application of the disordered local moment method to the low-spin phase of Pd3 Fe would begin to remedy the failings of density functional theory in this regime. Such calculations would be capable of testing the ideas presented at the end of Chapter 4 concerning the electronic structure changes underlying pressure-induced Invar behavior in Pd3 Fe. The crucial test would be to calculate 100 the spontaneous volume magnetostriction in Pd3 Fe at pressures around the pressure-induced Invar anomaly. If the application of pressure results in a large spontaneous volume magnetostriction in Pd3 Fe on par with the value in classical Invar Fe3 Pd at ambient pressure, then our hypothesis will be proven and a single theoretical framework will unite classical and pressure-induced Invar behavior in Fe-Pd. The use of the disordered local moment method to model the high-pressure magnetic phase in Pd3 Fe might improve the prediction of the magnetic transition pressure. However, the disordered local moment calculations used in Invar studies to this point rely on the local spin-density approximation (LSDA) [32]. Evidence exists that both the local spin-density approximation and the generalized gradient approximation are insufficient to unambiguously predict the ground state magnetic ordering in materials, such as Pd3 Fe, with many closely competing magnetic states [40]. Not only do the LDA and GGA often disagree as to the ground state magnetic structure at a given volume, but even different parameterizations within the GGA have been shown to disagree on magnetic ground states [40]. A solution to this problem would be to move to a many-body technique, which does not rely on the relatively simple approximations for exchange and correlation that the LDA and GGA employ (outlined in Chapter 3). A particularly promising direction in this regard is dynamical mean field theory (DMFT) [166, 167], which has already been applied successfully to Fe [169, 170, 171] and Fe-Ni [168] alloys. Due to its mean field-like nature [172], the disordered local moment method is also not capable of addressing the existence of local spin correlations in the high-pressure phase of Pd3 Fe suggested by the work presented in Chapter 5. More advanced ab initio spin dynamics simulations [173] might begin to address these challenging issues. The identification and characterization of new systems that are transformed into Invar materials only at high pressure will provide invaluable insight into the problems that remain outstanding around pressure-induced Invar behavior in Pd3Fe. With Fe-Ni and Fe-Pd already known to exhibit pressure-induced Invar behavior, the obvious next system to look into is Fe-Pt. Beyond this, all 101 known Invar systems should be examined for pressure-induced Invar behavior at compositions away from their classical Invar compositions. Additionally, the prospect of discovering pressure-induced Invar behavior in a system that does not display Invar behavior at any composition at ambient pressure would be particularly exciting. The efficiency with which density functional theory calculates electronic structures allows for computational screening of many alloy systems in much less time than it would take to carry out the high-pressure experimental characterization. Using the insights gained from the work presented in this thesis concerning the electronic structure correlates of pressure-induced Invar behavior, many alloy systems could be screened quickly using density functional theory methods, and the most promising selected for further experimental study. The ultimate value of such work would be a deeper understanding of Invar behavior in general, potentially leading to powerful new insights concerning itinerant electron magnetism and technological applications. We have taken a few steps into this new world of pressure-induced Invar behavior and found the terrain surprisingly rich. After more than a hundred years of searching, perhaps the insights gained by studying the pressure-induced Invar effect will finally cut the Gordian Knot of Invar complexity that has baffled so many for so long. With ever-advancing experimental and theoretical techniques available, and an ever-growing body of Invar knowledge to draw upon, may the Invar researchers of this new century accomplish that which those of the last century failed to accomplish, a complete microscopic understanding of the Invar effect. 102 Appendix A Flux, Brightness, and Brilliance In this appendix we build up to the definition of brilliance, following [52]. To begin, we define the total flux as, T otal F lux = P hotons . second (A.1) This definition includes photons of all energies, and so is of interest for measurements that use white radiation. To characterize monochromatic radiation we define the spectral flux, which accounts for the spectral distribution of radiation. The spectral flux is the flux normalized to a relative spectral bandwidth of ∆E/E = 10−3 [52], Spectral F lux = P hotons/s , 0.1%bandwidth (A.2) with s standing for second. Normalizing to a solid angle of 1 mrad2 we obtain the definition of brightness, Brightness = P hotons/s mrad2 0.1%bandwidth . (A.3) The definition of the Brilliance is arrived at by further accounting for the size of the x-ray source. Normalizing to a source area of 1 mm2 we have, 103 Brilliance = P hotons/s . mrad2 mm2 0.1%bandwidth (A.4) The Brilliance of a synchrotron source characterizes the photons per second available in a defined energy bandwidth, normalized to the solid angle of the radiation cone and the source area. 104 Appendix B Ab Initio Temperature Dependent Enthalpy of Reaction in Li-Mg-Si System Typically, density functional methods are considered zero-temperature, ground state techniques. The small-displacement method allows us to extend 0 K DFT results to finite temperature. The smalldisplacement method assumes the harmonic approximation, consequently the finite temperature thermodynamic results available from the small-displacement method also assume the harmonic (or quasi-harmonic) approximation. To see how the small-displacement method extends DFT results to finite temperature we will look at a specific example; the enthalpy of reaction (∆H) for the dehydrogenation reaction, MgH2 + 2LiH + Si → Li2MgSi + 2H2 . (B.1) This reaction is of interest for the storage of hydrogen for mobile fuel cell applications. A viable hydrogen storage material for automotive applications must meet stringent requirements for hydrogen storage capacity and weight. LiH and MgH2 , with 12.7 and 7.67 wt.% H2 storage capacity, respectively, have long been studied as potential hydrogen storage materials for mobile applications. The problem researchers encounter with light metal hydrides, such as LiH and MgH2 , is that hydrogen is bound too strongly in these materials to be released at temperatures conveniently reached using Proton Exchange Membrane fuel cell technology. An ideal hydrogen storage material for use with a 105 PEM fuel cell must show significant dehydrogenation below 100◦ C. To reach an equilibrium H2 pressure of 1 bar, MgH2 and LiH must be heated to 275◦C and 910◦C, respectively. LiH and MgH2 are too stable for practical mobile applications. Wiswall and Reilly [174, 175] pioneered the idea of destabilizing metal hydrides by adding constituents to the dehydrogenation reaction that stabilize the dehydrogenated state. This has the advantage lowering the dehydrogenation temperature, but comes at a cost of increased weight and so decreased H2 weight percent. Vajo et al. [176] destabilized the dehydrogenation reactions in LiH and MgH2 using Si to form stable Si-containing alloys in the dehydrogenated states. Other examples of increased H2 pressure at reasonable temperatures through destabilization have been studied, but to date none are fully satisfactory for mobile applications. Reaction B.1 has 5.73 wt.% H2 storage capacity, making it a viable hydrogen storage material if it releases significant quantities of hydrogen at sufficiently low temperatures. A proven estimate for the temperature of dehydrogenation is given by the relation [177], T = ∆H , ∆S (B.2) here T is the temperature of dehydrogenation, and ∆H and ∆S are the enthalpy and entropy changes caused by the dehydrogenation reaction. Alapati et al. [177] use [178] to estimate that for metal hydride reactions, 100 J J ≤ ∆S ≤ 130 . K K (B.3) This means that to obtain a dehydrogenation temperature between 300 and 600 K we need, 30 kJ mol−1 H2 ≤ ∆H ≤ 78 kJ mol−1 H2 . (B.4) DFT calculations allow us to estimate ∆H before experimental effort is expended on studying a system. Having a favorable ∆H is a necessary, but not sufficient, condition for any potentially viable hydrogen storage material. Even with a favorable ∆H a reaction can still be kinetically 106 limited, for example. In general, ∆H = X products H − X H. (B.5) reactants The method of calculating ∆H from first principles differs greatly for the solid crystals and the gaseous H2. We take the methods one at a time, beginning with our treatment of the solid phase reactants. So we begin with, ∆H(T ) = ∆Hs(T ) + ∆HH2 (T ), (B.6) with s denoting solids (M gH2, LiH, Si, Li2 M gSi). Now we can write, Hs(T ) = U (T ) + P V, (B.7) with U (T ) the temperature dependent internal energy of the solids, P the pressure and V the volume of the solids. At the temperatures of interest, ∆P V is negligible for the solids compared to the contributions of the other terms of interest; consequently ∆P V for the solid phases will be neglected throughout our analysis (following the lead of [177]). We rewrite Eq. B.7, neglecting the PV term and breaking the internal energy term into two parts, Hs(T ) = U0 + Uvib(T ). (B.8) U0 is the energy of the electronic system at 0 K as calculated by density functional theory. Uvib (T ) is the vibrational internal energy of the solids, originating from the vibrations of the atoms about their equilibrium lattice positions. To calculate U0 density functional theory calculations were carried out in the generalized gradient approximation [179] for the electronic exchange and correlation potential using the projector 107 augmented-wave (PAW) method [133, 104] as implemented in the Vienna Ab-initio Simulation Package (VASP) [90]. Unit cells for MgH2 (P42 /mnm, 6 atoms), LiH (Fm-3m, 8 atoms), Si (Fd-3m, 8 atoms), and Li2MgSi (P-43m, 16 atoms) were constructed. An energy cutoff between 250 and 370 eV was used together with an 8 × 8 × 8 Monkhorst–Pack grid [134]. Convergence with respect to both basis set cutoff energy and number of k-points was tested and found to be better than 2 meV/atom. To determine the ground state geometry for each material the forces on the atoms and the stresses on the unit cell were minimized by allowing the cell shape, cell volume and atomic positions of all atoms within the unit cell to relax simultaneously according to the conjugate gradient algorithm until the residual Hellmann–Feynman forces on all of the atoms were smaller than 0.1 meV/Å. An energy cutoff 1.3 times larger than the default energy cutoff was used during the relaxations to facilitate accurate calculation of the diagonal components of the stress tensor [92]. This resulted in the geometry and lattice parameters shown in Fig. B.1 and Table B.1. Moriga et al., proposed two possible structures for Li2MgSi [180], one with P-43m (No. 215) symmetry and one with Fm-3m (No. 225) symmetry. Ancillary calculations showed that the P-43m structure was the low energy phase at ambient pressure. Table B.1: Parameters obtained from DFT (GGA-PBE)calculations. Experimental numbers taken from [181] and references therein. E (eV/atom) a0 (Å) (DFT) a0 (Å) (Expt.) MgH2 -2.951 4.465 (c=2.994) 4.517 (c=3.021) LiH -3.050 3.981 4.085 Si -5.422 5.467 5.430 Li2 MgSi -2.932 6.373 6.415 [180] H2 -3.381 N/A N/A In the harmonic approximation the vibrational internal energy is calculated as [113], 1 Uvib(T ) = r 2 Z∞ ~ω coth  ~ω 2kB T  g(ω) dω, (B.9) 0 with Uvib (T ) the vibrational internal energy of a primitive unit cell, r the number of degrees of freedom in the unit cell, ω the vibrational frequency, and g(ω) the vibrational density of states. To obtain g(ω), supercells representing MgH2 (48 atom supercell), LiH (64 atom supercell), Si 108 Figure B.1: Crystal structures used for first-principles calculations on Li2MgSi (upper left), Si (upper right), MgH2 (lower left), and LiH (lower right). (64 atom supercell), and Li2MgSi (128 atom supercell) were constructed. In the supercells, one at a time, each symmetry inequivalent atom was displaced from its equilibrium position in each symmetry inequivalent direction. This entailed 22 total displacements. The Hellmann-Feynman forces induced on all other atoms from each displacement were determined from a total energy calculation performed on each displaced supercell using VASP. The resultant Hellmann-Feynman forces were used to generate the dynamical matrix using the direct method [110] as implemented in the PHONON program [113]. The diagonalization of the dynamical matrix at a series of points in reciprocal space gave the phonon dispersion curves and the phonon density of states. The resulting vibrational (total and partial) densities of states are shown in Fig. B.2. Using the calculated g(ω), 109 the integration in Eq. B.9 was carried out to obtain the temperature dependent vibrational internal energy for the reaction (shown in Fig. B.4). Figure B.2: Calculated partial and total phonon density of states for (a)MgH2 , (b)LiH, (c)Si, and (d)Li2MgSi. To treat the H2 gas we write [182, 183], HH2 (T ) = U0,H2 + UH2 (T ) + (P V )H2 = U0,H2 + Uvib(T ) + Utrans (T ) + Urot (T ) + (P V )H2 ! − hν 3 NA hν NA hνe kB T 2 + RT + RT + (P V )H2 = U0,H2 + + − khνT 2 2 2 1−e B hν ! − 7 NA hν NA hνe kB T + RT = U0,H2 + + hν − 2 2 1 − e kB T (B.10) with U0,H2 electronic energy of the H2 molecule calculated in DFT, NA Avogadro’s constant, h 110 Planck’s constant, ν the H2 vibrational frequency, kB the Boltzmann constant. The translational ((3/2)RT), rotational (RT), and pressure-volume (PV = RT, treating H2 as an ideal gas) contributions from the H2 gas are all accounted for in Eq. B.10. To find the electronic energy of the hydrogen molecule, H2 is placed in the center of a large box (12 Å×12 Å×12 Å). The size of the box is chosen so that images of the H2 molecule that result from the periodic boundary conditions inherent in the DFT code do not interact significantly with one another and distort our desired results. The H2 molecule was allowed to relax along its axis according to the conjugate-gradient algorithm until the residual Hellmann-Feynman forces on the atoms were smaller than 0.3 meV/Å, coming to an equilibrium separation of 0.751 Å. Figure B.3: Energy-separation data generated for the H atoms in the H2 molecule. The H2 vibrational frequency was obtained by varying the distance between the H atoms in the H2 molecule and performing a series of fixed-point energy calculations. The second derivative with respect to atomic separation of this energy-separation data (shown in Fig. B.3) was taken, giving the force constant for the H-H bond. Inserting this into, 111 1 ν= 2π r K m (B.11) with K the force constant and m the mass of a hydrogen atom gave the zero-point vibrational frequency of the H2 molecule. As shown in Fig. B.4, the enthalpy of reaction for, MgH2 + 2LiH + Si → Li2MgSi + 2H2 , (B.12) is 49.15 kJ/mol H2 at 0 K, 57.91 kJ/mol H2 at 300 K, and 59.68 kJ/mol H2 at 600 K. These results indicate that this dehydrogenation reaction is a viable candidate for use as a hydrogen storage material for mobile applications, i.e., it is likely to release significant amounts of hydrogen fuel at optimal temperatures for use with a PEM fuel cell. Experimental investigation of this reaction is underway at the present time. This example demonstrates how standard DFT results can be extrapolated to finite temperatures using the small-displacement method and standard thermodynmaics. 112 Figure B.4: (Upper) Calculated temperature dependent enthalpy of reaction for the reaction MgH2 + 2LiH + Si → Li2MgSi + 2H2 . (Lower) Calculated temperature dependence of the vibrational internal energy (Uvib(T)) for the reaction, the translational, rotational, and pressure-volume contributions from the H2 gas, as well as the 0 K electronic total energy of the reaction (U0 ). 113 Bibliography [1] L. Dubrovinsky, N. Dubrovinskaia, I. A. Abrikosov, M. Vennström, F. Westman, S. Carlson, M. van Schilfgaarde, and B. Johansson , Phys. Rev. Lett. 86, 4851 (2001). [2] A. Kussmann and K. Jessen, Phys. Soc. Japan Suppl. 17 B 1, 136 (1962). [3] R. J. Hemley and N. W. Ashcroft, Physics Today 51, 26 (1998). [4] C. E. Guillaume, C. R. Acad. Sci. 125, 235 (1897). [5] E. F. Wassermann, in Ferromagnetic Materials, edited by K. H. J. Buschow and E. P. Wohlfarth (North- Holland, Amsterdam, 1990), Vol. 5, p. 237. [6] R. J. Weiss, Proc. Phys. Soc. 82, 281 (1963). [7] V. L. Moruzzi, Phys. Rev. B 41, 6939 (1990). [8] P. Entel, E. Hoffmann, P. Mohn, K. Schwarz, and V. L. Moruzzi, Phys. Rev. B 47, 8706 (1993). [9] I. A. Abrikosov, O. Eriksson, P. Söderlind, H. L. Skriver, and B. Johansson, Phys. Rev. B 51, 1058 (1995). [10] M. Podgórny, Phys. Rev. B 46, 6293 (1992). [11] R. Hayn and V. Drchal, Phys. Rev. B 58, 4341 (1998). [12] C. A. Kuhnen et al., Phys. Rev. B 46, 8915 (1992). [13] P. Mohn, E. Supanetz, and K. Schwarz, Aust. J. Phys. 46, 651 (1993). [14] M. M. Abd-Elmedguid and H. Micklitz, Phys. Rev. B 40, 7395 (1989). 114 [15] S. Odin, F. Baudelet, J. P. Itié, A. Polian, S. Pizzini, A. Fontaine, C. Giorgetti, E. Prtyge, and J. P. Kappler, J. Appl. Phys. 83, 7291 (1998). [16] J. P. Rueff, A. Shukla, A. Kaprolat, M. Krisch, M. Lorenzen, F. Sette, and R. Verbeni, Phys. Rev. B 63, 132409 (2001). [17] G. Oomi and N. Mori, J. Phys. Soc. Japan 50, 291 (1981). [18] M. Matsushita, Y. Nakamoto, E. Xuzuki, Y. Miyoshi, H. Inoue, S. Endo, T. Kikegawa, and F. Ono, J. Magn. Magn. Mater. 284, 403 (2004). [19] F. Decremps and L. Nataf, Phys. Rev. Lett. 92, 157204 (2004). [20] L. Nataf, F. Decremps, M. Gauthier, and B. Canny, Phys. Rev. B 74, 184422 (2006). [21] S. Khmelevskyi, I. Turek, and P. Mohn, Phys. Rev. Lett. 91, 037201 (2003). [22] M. van Schilfgaarde, I. A. Abrikosov, and B. Johansson, Nature (London) 400, 46 (1999). [23] Y. Wang, G. M. Stocks, D. M. C. Nicholson. W. A. Shelton, V. P. Antropov, and B. N. Harmon, J. Appl. Phys. 81, 3873 (1997). [24] M. Uhl, L. M. Sandratskii, and J. Kübler, Phys. Rev. B 50, 291 (1994). [25] L. Mañosa, et al., Phys. Rev. B 45, 2224 (1992). [26] N. Cowlam and A. R. Wildes, J. Phys.: Condens. Matter 15, 521 (2003). [27] A. V. Ruban, S. Khmelevskyi, P. Mohn, and B. Johansson, Phys. Rev. B 76, 014420 (2007). [28] M. Matsushita, S. Endo, K. Miura, and F. Ono, J. Magn. Magn. Mater. 265, 352 (2003). [29] M. Matsushita, S. Endo, K. Miura, and F. Ono, J. Magn. Magn. Mater. 269, 393 (2004). [30] M. Matsushita, Y. Miyoshi, S. Endo, and F. Ono, Phys. Rev. B 72, 214404 (2005). [31] S. Khmelevskyi and P. Mohn, Phys. Rev. B 68, 214412 (2003). [32] S. Khmelevskyi and P. Mohn, Phys. Rev. B 69, 140404 (2004). 115 [33] S. Khmelevskyi, et al., Phys. Rev. B 72, 064510 (2005). [34] B. L. Gyorffy, et al., J. Phys. F: Met. Phys. 15, 1337 (1985); F. J. Pinski, et al., Phys. Rev. Lett. 56, 2096 (1986). [35] V. Crisan, P. Entel, H. Ebert, H. Akai, D. D. Johnson, J. B. Staunton, Phys. Rev. B 66, 014416 (2002). [36] E. F. Wassermann and M. Acet, in Magnetism and Structure in Functional Materials edited by A. Planes, et al. (Springer-Verlag, Berlin, 2005), p. 177. [37] M. Acet, E. F. Wassermann, K. Andersen, A. Murani, O. Schärpff, Europhys. Lett. 40, 93 (1997). [38] J. Kaspar and D. R. Salahub, Phys. Rev. Lett. 47, 54 (1981). [39] E. F. Wassermann, in Concise Encyclopedia of Magnetic and Superconducting Materials, edited by K. H. J. Buschow (Elsevier, Oxford, 2005), p. 336. [40] I. A. Abrikosov, A. E. Kissavos, F. Liot, B. Alling, S. I. Simak, O. Peil, and A. V. Ruban, Phys. Rev. B 76, 014434 (2007). [41] R. F. Sabirtyanov, S. K. Bose, and O. N. Mryasov, Phys. Rev. B 51, 8958 (1995). [42] P. Entel, E. Hoffmann, M. Clossen, K. Kadau, M. Schröter, R. Meyer, H. C. Herper, and M. S. Yang, in The Invar Effect: A Centennial Symposiu, edited by J. Wittenauer (Metallurgical Society, Warrendale, 1996), p. 87-104. [43] P. Entel, H. C. Herper, E. Hoffmann, G. Nepecks, E. F. wassermann, M. Acet, V. Crisan, and H. Akai, Phil. Mag. B 80, 141 (2000). [44] P. J. Brown, K. U. Numann, and K. R. A. Ziebeck, J. Phys.: Condens. Matter 13, 1563 (2001). [45] P. J. Brown, T. Kanomata, M. Matsumoto, K. U. Neumann, and K. R. A Ziebeck, J. Magn. Magn. Mater. 242, 781 (2002). 116 [46] P. J. Brown, T. Chaterji, J. Kaestner, and K. R. A. Ziebeck, J. Phys.: Condens. Mat. 16, 5249 (2004). [47] L. Merrill and W. A. Bassett, Rev. Sci. Instum. 45, 290 (1974). [48] E. Sterer, M. P. Pasternak, and R. D. Taylor, Rev. Sci. Instr. 61, 1117 (1990). [49] R. A. Formman, G. J. Piermarini, J. D. Barnett, and S. Block, Science 176, 284 (1972). [50] J. P. Poirier, Introduction to the Physics of the Earth’s Interior (Cambridge University Press, Cambridge, 2000), Second Edition. [51] F. D. Murnaghan, Proc. Natl. Acad. Sci. 30, 244 (1944). [52] R. Röhlsberger, Nuclear Condensed Matter Physics with Synchrotron Radiation (SpringerVerlag, Berlin, 2004). [53] K. Balewski, W. Brefeld, W. Decking, W. Drube, H. Franz, R. Gehrke, P. Gürtler, U. Hahn, J. Pflüger, H. Schulte-Schrepping, M. Tischer, and E. Weckert, DESY-HASYLAB Internal Report (2002). [54] Argonne National Laboratory, Lightsources Image Bank, http://www.lightsources.org/imagebank . [55] B. Fultz and J. M. Howe, Transmission Electron Microscopy and Diffractometry of Materials (Springer-Verlag, Berlin, 2002), Second Edition. [56] W. L. Bragg, Proc. Camb. Phil. Soc. 17, 43 (1913). [57] H. M. Rietveld, Acta Crystallogr. 22, 151 (1967). [58] R. L. Mössbauer, Z. Physik 151, 124 (1958). [59] B. Fultz, in Methods in Materials Research: A Current Protocols Publication edited by B. Fultz, et al. (John Wiley, New York, 2000), Unit 9.c.1. [60] A. B. Papandrew, Ph.D. diss., California Insititute of Technology, 2006. 117 [61] Argonne National Laboratory, http://www.aps.anl.gov/Facility/StorageRingParameters. [62] S. L. Ruby: J. Phys. Colloques 35, C6-209 (1974). [63] E. Gerdau, RR. Rüffer, H. Winkler, W. Tolksdorf, C. P. Klages, and J. P. Hannon: Phys. Rev. Lett. 54, 835 (1985). [64] W. Sturhahn, T. S. Toellner, E. E. Alp, X. Zhang, M. Ando, Y. Yodo, S. Kikuta, M. Seto, and C. W. Kimball, Phys. Rev. Lett. 74, 3832 (1995). [65] Y. Seto, Y. Yoda, S. Kikuta, X. W. Zhang, and M. Ando, Phys. Rev. Lett. 74, 3828 (1995). [66] A. I. Chumakov. R. Rüffer, H. Grünsteudel, H. F. Grünsteudel, G. Grübel, J. Metge, O. Leupold, and H. A. Goodwin, Europhys. Lett. 30, 427 (1995). [67] L. Zhang, J. Stanek, S. S. Hafner, H. Ahsbahs, H. F. Grünsteudel, J. Metge, and R. Rüffer, American Mineralogist 84, 447 (1999). [68] W. Sturhahn, Hyperfine Interact. 125, 149 (2000). [69] M. Y. Hu, W. Sturhahn, T. S. Toellner, P. M. Hession, J. P. Sutter, E. E. Alp, Nucl. Instrum. Meth. A 428, 551 (1999). [70] Y. L. Chen, D. P. Yang, Mössbauer Effect in Lattice Dynamics (WILEY-VCH, Weinheim, 2007). [71] D. J. Griffiths, Introduction to Quantum Mechanics (Prentice Hall, Upper Saddle River, 1995). [72] R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, Cambridge, 2004). [73] C. M. Bakewell, Source Book in Ancient Philosophy (Charles Scribner’s Sons, New York, 1907). [74] W. Kohn, Rev. Mod. Phys. 71, 1253 (1999). [75] N. W. Ashcroft and N. D. Mermin, Solid State Physics (Holt, Rinehart and Winston, New York, 1976). [76] M. Born and J. R. Oppenheimer, Ann. Physik 84, 457 (1927). 118 [77] F. Seitz, Modern Theory of Solids (McGraw-Hill, New York, 1965). [78] M. Born and K. Huang, Dynamical Theory of Crystal Lattices (Clarendon, Oxford, 1954). [79] J. C. Slater, Quantum Theory of Molecules and Solids (McGraw-Hill, New York, 1963). [80] J. Ziman, Principles of the Theory of Solids (Cambridge University Press, London, 1965). [81] G. Venkataraman, L. A. Feldkamp, and V. C. Sahni, Dynamics of Perfect Crystals (MIT Press, Cambridge, 1975). [82] L. H. Thomas, Proc. Camb. Phil. Soc. 23, 542 (1927). [83] E. Fermi, Atti. Accad. Nazl. Lincei 6, 602 (1927). [84] P. Hohenberg and W. Kohn, Phys. Rev. B 136, 864 (1964). [85] W. Koch and M. C. Holthausen, A Chemist’s Guide to Density Funcational Theory (WileyVCH, Einheim, 2000). [86] M. D. Towler, Phys. Stat. Sol. (b) 243, 2573 (2006). [87] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). [88] S. Redner, Phys. Today 58, 49 (2005). [89] E. Schrödinger, Phys. Rev. 28, 1049 (1926). [90] G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996); Phys. Rev. B. 54, 11169 (1996). [91] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in FORTRAN (Cambridge University Press, Cambridge, 1992). [92] G. Kresse, M. Marsman, and J. Furthmüller, VASP the Guide (Universität Wien, Wien, 2009). [93] E. R. Davidson, in Numerical Algorithms in Chemistry: Algebraic Methods, edited by C. Moler and I. Shavitt (Univ. of California, Berkeley, 1978), p. 15. 119 [94] E. R. Davidson, in Methods in Computational Molecular Physics, edited by G. H. F. Diercksen and S. Wilson, (Plenum, New York, 1983), p. 95. [95] E. R. Davidson, J. Comput. Phys. 17, 87 (1975). [96] B. Liu, in Numerical Algorithms in Chemistry: Algebraic Methods, edited by C. Moler and I. Shavitt (Univ. of California, Berkeley, 1978),p. 49. [97] M. P. Teter, M. C. Payne, and D. C. Allan, Phys. Rev. B 40, 12255 (1989). [98] D. M. Bylander, L. Kleinman, and S. Lee, Phys. Rev. B 42, 1394 (1990). [99] P. Pulay, Chem. Phys. Lett. 73, 393 (1980). [100] D. M. Wood and A. Zunger, J. Phys. A: Math. Gen. 18, 1343 (1985). [101] D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980); J.P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). [102] J. Thijssen, Computational Physics (Cambridge University Press, New York, 2007), Second Edition. [103] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994). [104] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). [105] R. D. King-Smith and R. J. Needs, J. Phys.: Condens. Matter 2, 3431 (1990). [106] K. Karch, P. Pavone, W. Windl, O. Schütt, and D. Strauch, Phys. Rev. B 50, 17054 (1994). [107] O. Schütt, P. Pavone, W. Windl, K. Karch, and D. Strauch, Phys. Rev. B 50, 3746 (1994). [108] K. Kunc and R. M. Martin, Phys. Rev. Lett. 48, 406 (1982). [109] W. Frank, C. Elsässer, and M. Fähnle, Phys. Rev. Lett. 74, 1791 (1995). [110] G. Kresse, J. Furthermüller, and J. Hafner, Europhys. Lett. 32, 729 (1995); K. Parlinski, Z. Q. Li, and Y. Kawazoe, Phys. Rev. Lett. 78, 4063 (1997); K. Parlinski, Z. Q. Li, and Y. Kawazoe, Phys. Rev. B 61, 272 (2000). 120 [111] H. Hellmann, Einführung in die Quantenchemie (Deuticke, Lepzig, 1937). [112] R. P. Feynman, Phys. Rev. 56, 340 (1939). [113] K. Parlinski, PHONON Manual, Version 4.24 (Polish Academy of Sciences, Cracow, 2005). [114] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes (Cambridge: Uhiversity Press, Cambridge, 1992), p. 670. [115] R. J. Jastrow, Phys. Rev. 98, 1479 (1955). [116] W. M. C. Foulkes, M. Nekovee, R. L. Gaudoin, M. L. Stedman, R. J. Needs, R. Q. Hood, G. Rajagopal, M. D. Towler, P. R. C. Kent, Y. Lee, W. K. Leung, A. R. Porter, and S. J. Breuer, in High-Performance Computing, edited by R. J. Allan, M. F. Guest, A. D. Simpson, and D. S. Henty (Plenum, New York, 1999), p. 165. [117] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). [118] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001). [119] J. B. Anderson, J. Chem. Phys. 65, 4121 (1976). [120] R. J. Needs, M. D. Towler, N. D. Drummond, and P. Loppez Rios, CASINO version 2.0 User Manual (University of Cambridge, Cambridge, 2006). [121] N. D. Drummond and R. J. Needs, Phys. Rev. B 73, 024107 (2006). [122] X. Gonze, J.-M. Beuken, R. Caracas, F. Detraucx, M. Fuchs, G. -M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, Ph. Ghosez, J. -Y. Raty, and D. C. Allan, Computational Materials Science 25, 478 (2002). [123] D. D. Johnson, et al., in Physical Metallurgy of Controlled Expansion Invar-Type Alloys, edited by K. C. Russel and D. F. Smith (TMS, Warrendale, PA 1990), p. 3; I. A. Abrikosov, et al., Phys. Rev. B 76, 014434 (2007). [124] T. Sohmura, et al., Scr. Metall. 14, 855 (1980); R. Oshima, Scr. Metall. 19, 315 (1985). 121 [125] J. Crangle and W. R. Scott, J. Appl. Phys. 36, 921 (1965); G. G. Low, Adv. Phys. 18, 371 (1969). [126] Y. Tsunoda and R. Abe, Phys. Rev. B 55, 11507 (1997); Y. Tsunoda and R. Abe, Physica B 237-238, 458 (1997); R. Abe, et al., J. Phys.: Condens. Matter 1998, L79 (1998). [127] A. J. Smith, et al., J. Phys. F 7, 2411 (1977). [128] G. Longworth, Phys. Rev. 172, 572 (1968). [129] M. Fallot, Ann. Physique 10, 291 (1938). [130] S. J. Pickart and R. Nathans, J. Appl. Phys. 33, 1336 (1962); V. A. Tsurin and A. Z. Men’shikov, Phys. Met. Metall 45, 82 (1978); A. Z. Men’shikov and V. A. Tsurin, Phys. Met. Metall. 47, 68 (1979). [131] A. Jayaraman, Rev. Modern Physics 55, 65, (1983). [132] M. Hansen, Constitution of Binary Alloys, Metallurgy and Metallurgical Engineering Series (McGraw-Hill, New York, 1958), Second Edition. [133] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994); G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). [134] H.J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976). [135] M. Methfessel and A. T. Paxton, Phys. Rev. B 40, 3616 (1989). [136] C. A. Kuhnen, et al., Phys. Rev. B 35, 370 (1987); C. A. Kuhnen, et al., Phys. Rev. B 46, 8915 (1992); S. K. Bose, et al., J. Magn. Magn. Mat. 87, 97 (1990); T. Nautiyal, et al., Physica Scripta 46, 527 (1992), Y. S. Shi, et al., Phys. Rev. B 65, 172410 (2002). [137] J. W. Cable, et al., Phys. Rev. 138, A755 (1965). [138] Sergii Khmelevskyi and Peter Mohn, Phys. Rev. B 68, 214412 (2003). [139] P. Gorria, et al., Phys. Status Solidi (RRL) 3, 115 (2009). 122 [140] P. Mohn, Magnetism in the Solid State, Springer Series in Solid State Sciences Vol. 134 (Springer-Verlag, Berlin, Heidelberg, 2003). [141] S. Khmelevskyi, et al., Phys. Rev. Lett. 91, 037201 (2003); S. Khmelevskyi and P. Mohn, Phys. Rev. B 69, 140404 (2004); S. Khmelevskyi, et al., Phys. Rev. B 72, 064510 (2005). [142] K. Tanaka and K. Morioka, Phil. Mag. 83, 1797 (2003). [143] T. Sohmura et al., Scr. Metall. 14, 855 (1980); R. Oshima, Scr. Metall. 19, 315 (1985). [144] A. Kussmann and K. Jessen, J. Phys. Soc. Japan Suppl. 17B1, 136 (1962). [145] W. G. Stirling, R. A. Cowley, and M. W. Stringfellow, J. Phys. F: Metal Phys. 2, 421 (1972). [146] A. F. Yue, A. Papandrew, P. D. Bogdanoff, I. Halevy, J. G.-W. Lin, B. Fultz, W. Sturhahn, E. E. Alp, and T. S. Toellner, Hyperfine Interactions 141/142, 249 (2002). [147] S. Ghosh, J. Phys.: Condens. Matter 20, 275208 (2008). [148] S. Baroni, S. De Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001). [149] J. M. Leger, C. Susse-Loriers, and B. Vodar, Phys. Rev. B 6, 4250 (1972); G. Hausch, Phys. Status Solidi 16, 371 (1973); G. Oomi and N. Mōri, J. Phys. Soc. Japan 50, 1043 (1981); J. Phys. Soc. Japan 50, 2917 (1981); J. Phys. Soc. Japan 50, 2924 (1981); M.M. Abd-Elmeguid, B. Schleede, and H. Micklitz, J. Magn. Mater. 72, 253 (1988); Phys. Rev. B 40, 7395 (1989). [150] S. Odin, F. Baudelet, J. P. Itié, A. Polian, S. Pizzini, A. Fontaine, Ch. Giorgetti, E. Dartyge, J. P. Kappler, J. Appl. Phys. 83, 7291 (1998); S. Odin, F. Baudelet, Ch. Giorgetti, E. Dartyge, J. P. Itié, A. Polian, J. C. Chervin, S. Pizzini, A. Fontaine, and J. P. Kappler, Europhys. Lett. 47, 378 (1999). [151] J. P. Rueff, A. Shukla, A. Kaprolat, M. Krisch, M. Lorenzen, F. Sette, and R. Verbeni, Phys. Rev. B 63, 1324091 (2001); M. Matsushita, T. Nishimura, S. Endo, M. Ishizuka, K. Kindo, and F. Ono, J. Phys.: Condens. Matter 14, 10753 (2002); M. Matsushita, S. Endo, K. Miura, and 123 F. Ono, J. Magn. Mater. 260, 371 (2003); M. Matsushita, S. Endo, K. Miura, F. Ono, J. Magn. Mater. 269, 393 (2004). [152] M. Matsushita, Y. Nakamoto, E. Suzuki, Y. Miyoshi, H. Inoue, S. endo, T. Kikegawa, F. Ono, J. Magn. Mater. 284, 403 (2004); M. Matsushita, Y. Miyoshi, S. Endo, and F. Ono, Phys. Rev. B 72, 2144041 (2005); F. Decremps and L. Nataf, Phys. Rev. Lett. 92, 1572041 (2004); L. Nataf, F. Decremps, M. Gauthier, and B. Canny, Phys. Rev. B 74, 1844221 (2006). [153] M. L. Winterrose, M. S. Lucas, A. F. Yue, I. Halevy, L. Mauger, J. A. Muñoz, Jingzhu Hu, M. Lerche, and B. Fultz, Phys. Rev. Lett. 102, 237202 (2009). [154] H. K. Mao, J. Xu, V. V. Struzhkin, J. Shu, R. J. Hemley, W. Sturhahn, M. Y. Hu, E. E. Alp, L. Vocadlo, D. Alf, G. D. Price, M. J. Gillan, M. Schwoerer-Bhning, D. Husermann, P. Eng, G. Shen, H. Giefers, R. Lbbers, G. Wortmann, Science 292, 914 (2001). [155] M. Born and K. Huang, Dynamical Theory of Crystal Lattices (Clarendon, Oxford, 1988). [156] A. Maradudin, E. Montroll, G. Weiss, and I. Ipatova, in Solid State Physics, 2nd ed., edited by H. Ehrenreich, F. Seitz, and D. Turnbull (Academic, New York, 1971), Suppl. 3. [157] K. V. Price, R. M. Storn, and J. A. Lampienen, Differential Evolution (Birkhäuser, 2005). [158] M. McKerns, P. Hung, and M. Aivazis, “Mystic: a simple model-independent inversion framework”, http://dev.danse.us/trac/mystic (unpublished). [159] M. G. Kresch, Ph.D. diss., California Institute of Technology, 2009. [160] R. Abe, Y. Tsunoda, M. Nishi, Nd K. Kakurai, J. Phys.: Condens. Matter 10, L79 (1998). [161] Y. Tsunoda and R. Abe, Phys. Rev. B 55, 11507 (1997). [162] Y. Tsunoda and R. Abe, Physica B 237-238, 458 (1997). [163] K. Fushimi and Y. Tsunoda, J. Phys. Soc. Jpn. 70, 1364 (2001). [164] M. Hirano, Y. Tsunoda, Phys. Rev. B 59, 13835 (1999). 124 [165] Y. Tsunoda, N. Hiruma, J. L. Robertson, and J. W. Cable, Phys. Rev. B 56, 11051 (1997). [166] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996). [167] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006). [168] J. Minar, L. Chioncel, A. Perlov, H. Ebert, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 72, 045125 (2005). [169] A. I. Lichtenstein, M. I. Katsnelson, and G. Kotliar, Phys. Rev. Lett. 87, 067205 (2001). [170] L. Chioncel, L. Vitos, I. A. Abrikosov, J. Kollar, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 67, 235106 (2003). [171] A. Grechnev, I. Di Marco, M. I. Katsnelson, A. I. Lichtenstein, J. Wills, and O. Eriksson, Phys. Rev. B 76, 035107 (2007). [172] S. Khmelevskyi, A. V. Ruban, and P. Mohn, J. Phys.: Condens. Matter 17, 7345 (2005). [173] V. P. Antropov, M. I. Katsnelson, M. van Schilfgaarde, and B. N. Harmon, Phys. Rev. Lett. 75, 729 (1995). [174] J. J. Reilly and R. H. Wisswall, Inorg. Chem. 6, 2220 (1967). [175] J. J. Reilly and R. H. Wiswall, Inorg. Chem. 7, 2254 (1968). [176] J. J. Vajo, F. Mertens, C. C. Ahn, R. C. Bowman, Jr. , and B. Fultz, J. Phys. Chem. B 108, 13977 (2004). [177] S. V. Alapati, J. K. Johnson, D. S. Sholl, Phys. Chem. Chem. Phys. 9, 1438 (2007). [178] A. Züttel, P. Wenger, S. Rentsch, P. Sudan, P. Mauron, and C. Emmenegger, J. Power Sources 118, 1 (2003). [179] J. P. Perdew, K. Bourke, M. Ernzerhof, Phys. Rev. Lett. 77 3865 (1996); Phys. Rev. Lett. 78, 1396 (1997). 125 [180] T. Moriga, K. Watanabe, D. Tsuji, S. Massaki, and I. Nakabayashi, Journal of Solid State Chemistry 153, 386 (2000). [181] S. V. Alapati, J. K. Johnson, and D. S. Sholl, J. Phys. Chem. B 110, 8769 (2006); J. Phys. Chem. C 111, 1584 (2007), Phys. Chem. Chem. Phys. 9, 1438 (2007); Phys. Rev. B 76, 104108 (2007). [182] D. A. McQuarrie and J. D. Simon, Molecular Thermodynamics (University Science Books, Sausalito, 1999). [183] D. A. McQuarrie Statistical Mechanics (University Science Books, Sausalito, 2000).