account of metallic nanocontacts and scalar-relativistic effects

Get Complete Project Material File(s) Now! »

OVERVIEW OF STANDARD METHODS

For convenience, the main modelling techniques used in this work are described in this chapter. A basic understanding of these techniques is essential in order to appreciate the developments and results that are presented in subsequent chapters.
Section 3.1 gives an overview of classical molecular dynamics (CMD), which is a well-established simulation technique dating back to the 1950s [146]. More specifically, Sec 3.1.1 describes the embedded-atom method (EAM) formalism, which is a model of reactive bonding between atoms in CMD simulations. It also briefly introduces the modified embedded-atom method (MEAM) formalism in which covalency is added to the isotropic bonding modelled by ordinary EAM potentials. Then, Sec 3.1.2 introduces an attempt to add collinear magnetism to CMD simulations within the EAM formalism.
Section 3.2 develops the more recent idea of spin-lattice dynamics (SLD), which was developed by Dudarev and co-workers in the late 2000s [48]. In the SLD formalism, there is an additional interaction between the atoms that depends on the relative noncollinear orientations of spins, which are treated semi-classically.
Section 3.3 deals with the fully quantum mechanical approach of density functional theory (DFT). Already formulated in the 1960s, by Kohn and Sham [147], DFT is today a widely-used simulation technique that has been considerably refined. Unlike classical molecular dynamics, which does not take into account electrons, DFT has the advantage of producing the full electronic band structure of the modelled materials, although this additional information does come at a much higher computational cost.
Finally, in Sec 3.4, the non-equilibrium Green’s function (NEGF) approach, as applied to DFT transport problems, is briefly described. The NEGF formalism is needed to treat out-of-equilibrium quantum mechanical systems, which usually have open boundaries and lack translation invariance. Like DFT, it was formulated in the 1960s, independently by Kadanoff and Bohm [148] and Keldysh [149], respectively.
Since none of the material that is reviewed in this chapter is fundamentally new, readers who are already familiar with the above-mentioned techniques may wish to skip certain sections, or simply proceed directly to the next chapter.

Classical molecular dynamics

CMD simulations are presently considered a standard method with well defined limitations, applications and validation tools. In the following sections the main aspects of this method are discussed.
The evolution of the atomistic structure of materials is widely modelled by means of CMD simulations, which makes systems involving millions of atoms now accessible [79]. In general, this approach involves solving Newton’s second law [150] for each of the particles in the system, as they evolve in time:
A distinct advantage of CMD simulations over other approaches, such as quantum mechanical molecular mechanics (QM/MM), is that the interaction potential V between the particles is treated semi-classically, and can usually be truncated, e.g., for electrically neutral systems such as metals, at an appropriately chosen cut-off separation between the atoms [146]. The value of the cut-off is usually fixed by fitting the interaction potential to ab-initio calculations and experimental properties of the system under consideration [146]. This greatly reduces the computational overhead required to model the system of particles, because the sum over interactions between all particles to calculate the total energy of the system, is limited by the cut-off.
In most applications, CMD models the properties of a macroscopic-sized sample of atomic and molecular systems. Unfortunately, the choice of boundaries (e.g., fixed or free or periodic) does not guarantee those properties are not affected by the size of the system chosen to represent the macroscopic sample. A standard illustration of the importance of simulation domain size in simulations of a three-dimensional system containing N particles, is that the particles at the surface of the domain number approximately N 1/3 [146]. Thus, in a simple cubic (SC) lattice with 1000 atoms, with roughly 49% of all the atoms located at the surface, free boundaries will constitute a very bad model of a macroscopic-sized sample [146]. The edge, or finite-size, effects that inevitably results from using free boundaries, are, however, greatly reduced in an SC lattice containing a million atoms, say, since only 6% of the atoms are now located at the surface of the simulation domain.
Periodic boundary conditions therefore offer an appropriate solution to the problem created by free boundaries, because they effectively mimic macroscopic structure. In practice, with periodic boundaries in force, the simulation domain is akin to a primitive cell within an infinite periodic lattice of identical cells (See Fig. 3.1) [146]. The red dotted line in Fig. 3.1 shows that it is important to ensure that periodic images of the atom in the central cell across the infinite simulation domain, do not overlap directly with the atom in the original cell. Otherwise, when a sum over energies is performed to calculate the total energy of the system, its value could become unrealistically large.
Periodic and shrink-wrapped boundaries are used in the CMD simulations in this work. However, overlap between periodic image atoms in contiguous cells or atoms on the opposite sides of central cells are sometimes provoked by the geometries of certain initial structures and actions applied to structures during simulations of nanocontact evolution. LAMMPS provides an easy way to avoid this problem by allowing for the use of shrink-wrapped boundary conditions. Unlike fixed boundary conditions, no atoms are lost when the simulation system crosses the boundary with shrink-wrapped boundary conditions in force. Instead the boundaries move with the simulation system when it exceeds the initial confines of the simulation domain.
Fig. 3.1: A two-dimensional simulation box containing three atoms (black square in the centre), which, in turn, has been embedded in an infinite two-dimensional lattice via periodic boundary conditions. The small black square encloses the minimum possible cell size that can be repeated in all directions. The red dotted line shows a periodic image from a contiguous cell, of the light green atom in the central cell. The black arrows represent pairwise interactions between atoms in contiguous smallest cells.6

The embedded-atom method potentials

There are many systems for which quantum mechanical effects are relatively unimportant, and for which CMD can give surprisingly accurate results. However, in cases where the classical approximation is valid, the success or failure of the method may still rely critically on the accuracy of the interaction potentials that are used in the model.
The class of potentials that can model the atomistic evolution of metallic nanocontacts require the ability to accurately reflect the reorganization and breaking of bonds between metal atoms. Potentials that can model bond making and breaking are referred to as reactive potentials [79]. They can vary in sophistication from the simple embedded-atom method (EAM) potential [31,78,79], used to model the metals in this thesis, to ReaxFF potentials [151], not applied here, which permit chemical reactions involving (originally) hydrocarbons.
In classical molecular dynamics within the EAM formalism, the potential energy of a system of N particles can be expressed as:
The EAM formalism of isotropic bonding between metal atoms provides a very realistic description of the structures adopted by metals, even in nanocontacts at first- and last-contact, as evidenced by the myriad previously reported results (see Refs. [11,17] and references therein). However, isotropic bonding remains a deficiency of the EAM model, and the so-called modified embedded-atom method (MEAM) potential [138] represents an attempt to correct this deficiency by including higher order moments of the electron “density” in the EAM formalism in Eq. (3.3). This approach lends directionality, or covalency, to the bonding between the atoms.
MEAM potentials nevertheless suffer from some of the same drawbacks as EAM potentials, in that they are typically fitted to low temperature properties of the materials they are intended to model, T = 0 K for properties obtained in ab initio calculations, [45,79]. Recently, a MEAM potential has been developed for Fe, Ni, Al and Cu that has not only been fitted to the melting point of these metals, but also to their near-melting point elastic constants. It thus very accurately reproduces the behaviour of the metals near their melting points [45], such as the interface between the solid and molten phases of a bulk metal. This is useful for the simulations performed in this work, because nanocontacts go through successive cycles of elastic and plastic deformation when they are subjected to cyclic loading. Of the CMD codes used in this work, the MEAM formalism is only available in the Large-scale atomic/molecular massively parallel simulator (LAMMPS), and will be used to model Ni and Fe nanocontacts, where the difference between directional vs. non-directional bonding may be important.

Magnetic potentials

CMD simulations can provide, for example, kinetic and potential energies, atomic stresses and forces, and atomistic (structural) details, but no electric or magnetic properties. However, in 2005, a physically reasonable ‘magnetic’ interatomic potential for CMD simulations of BCC iron was developed by Derlet and Dudarev [153]. It is an EAM-type potential for which ρi in Eq. (3.3), assigned to each atom i , determines the magnetic moment ζi of atom i , in a similar manner as functionals in DFT depend on the electron density at a given position in the system. We briefly consider the most important aspects here.
Derlet and Dudarev [153] used an approximation of the total energy per atom, Etot , as a function of its atomic magnetic moment ζ , based on an analogy with the simplest model of the second-order magnetic phase transition: the Ginzberg-Landau model (see Refs. [153,154] and references therein):
where α < 0 and β > 0 are fitting coefficients that depend on the local environment of each atom, and E0 is the energy of the non-magnetic phase [153,154]. Eq. (3.5) has the shape of a double-well potential with two energy minima. The roots of these minima are equal in magnitude but opposite in sign, and correspond to the equivalent “spin-up” and “spin-down” collinear states in an equilibrium ferromagnetic arrangement [153,154]:
The above model inspired Derlet and Dudarev to use an upside-down parabola to approximate the density of states D(E ) of the non-magnetic phase near the Fermi energy (centred at zero for convenience). Derlet and Dudarev parameterised D(E ) in terms of the d-orbital bandwidth W and two fitting parameters a > 0 and b > 0 that are independent of W :
The function F , chosen to be scalable so that D(E ) is independent of W , is a sum of parabolic and regular parts [153,154], the latter becoming negligible near the Fermi energy, EF≈0 , i.e., |R (0)|≪a . The total per-atom energy can also be written as [153]
It is well-known that the bandwidth W of the d-orbitals of Fe decreases as these orbitals become more localised, i.e., when the interatomic separation between the Fe atoms increases, and/or the coordination about every atom decreases. At the same time, a narrower bandwidth implies a higher density of states at the Fermi level D(E F) in 3d transition metals because the number of electrons in the d-bands must remain conserved. The famous Stoner criterion [155] predicts that if the product of the Stoner parameter I of the 3d metal and its density of states at the Fermi energy D(E F) exceeds unity, the material will be ferromagnetic. One can thus arrive at an expression for the EAM embedding functional F [ρi ] in Eq. (3.2) if one assumes the bandwidth W is proportional to the function ρi representing an effective “electron density”:
where A , B and ν are constants, and ρc is a critical effective “electron density” below which Fe is ferromagnetic (because the bandwidth is narrower then), while above ρc it is non-magnetic (much wider bandwidth).
To avoid the cusp at ρ≈ρc in Eq. (3.13), and thus ensure its first and second derivatives are continuous –in this way meeting the criteria for a second-order phase transition discussed earlier– the following final expression was adopted by Derlet and Dudarev [153,154]:
where C = 2.929μB and γ = 0.259 , obtained from a fit of Eq. (3.15) to DFT calculations of the magnetic moment as a function of the volume per atom in bulk iron [154]. The values of C and γ lead to an equilibrium per-atom magnetic moment of 2.154 μB for bulk BCC iron [154]. The agreement with the experimental value of 2.12 μB for the spin-only magnetic moment, or in DFT calculations in the local spin density approximation (LSDA), 2.15 μB , are both rather good [63].

Spin-lattice dynamics

In the last section, it was seen that the magnitude of the magnetic moment of an atom varies depending on its local environment in an assembly of ferromagnetic atoms. Lower coordination about a given ferromagnetic atom leads to greater confinement of its 3d electrons, which, in turn, enhances its magnetic moment (more unpaired electrons) because of greater Coulombic repulsion between the electrons.
This simple picture, which is derived from the Stoner formalism discussed in the previous section, is a model of long-range magnetism, finding its origin in the interplay between intra-atomic exchange and interatomic quantum hopping of valence electrons [48]. Both Ni and Fe exhibit long-range magnetism because of their delocalised 3d electrons. This fact makes these metals notoriously difficult to model at finite temperatures, even in quantum approaches [54,156–158].
Intuitively, however, the lower coordination of first neighbours (8 of them) about a given Fe atom in a bulk BCC crystal lattice, implies that there should be less overlap between its 3d orbitals than in Ni, which has 12 first neighbours in an FCC crystal [158]. Because the 3d electrons should, in the case of Fe, be relatively more localised than in Ni, a model of localised interacting magnetic moments, such as a generalised Heisenberg model of ferromagnetic exchange [48–50], may represent a reasonable description of the magnetism in Fe [54,157].
In fact, in a semi-empirical approach, the exchange energy of a weakly-inhomogeneous spin-polarisation density should always be modelled by a Heisenberg Hamiltonian (see Refs. [159–161] and references therein). Furthermore, at 4.2 K, the operating temperature in many STM and MCBJ experiments, longitudinal excitations of the magnetic moments in Fe, i.e., fluctuations in their magnitudes, are far less important than transverse excitations, i.e., precession of the magnetic moments about a local effective magnetic field produced by the moments on neighboring atoms [63,156].
It therefore becomes clear that the Derlet-Dudarev interatomic potential discussed in the previous section, which models collinear magnetism in Fe at 0 K in CMD simulations [48], cannot account for non-collinear spins or how these spins interact locally with each other [48], since they are modelled by an interatomic potential which allows only for spin “up” or spin “down”, and not any other orientation. In addition, the conservation of total angular momentum in ferromagnetic materials and, correspondingly, how energy and angular momentum are dissipated in spin currents, may lead to important technological ramifications in spintronics (See references in Ref. [48]). Conservation of total angular momentum in ferromagnetic materials requires non-collinear spin configurations and coupling between the spin and lattice degrees of freedom [55,56]. Therefore, at least for metals such as Fe, it is necessary to model the evolution of the three Cartesian spin degrees of freedom (Six , S iy , Siz ) of each atom i in socalled spin-lattice dynamics (SLD).
In this regard, a semi-classical exchange model of non-collinear magnetism was implemented for bulk Fe in 2008 by Ma et al. [48], and released in 2016 as a freely available spin-lattice dynamics code SPILADY [134]. The exchange parameters used in this model were procured by means of a standard approach followed in spin-polarised DFT, whereby ab-initio data are mapped onto classical Hamiltonians, such as a generalised Heisenberg exchange Hamiltonian [54,63]. The SLD model in Ref. [48] was developed out of interest in the effects lattice vibrations have on the stability of magnetism in FCC versus BCC iron [52,162]. To model such effects, Ma et al. added a generalised Heisenberg exchange term to the Derlet-Dudarev potential, ΗDD = EEAM in Eq. (3.2) [48]:
where J 0 = 904.90177 meV, and the cut-off radius r c = 3.75 Å was chosen to lie between 2nd and 3rd nearest neighbors in BCC Fe [48]. Θ(r c−r ij) is again the Heaviside step function.
The choice of the simple isotropic function in Eq. (3.17) was motivated by the small effects lattice vibrations have on magnons in BCC iron [52,54], and also by the fact that forces in a molecular dynamics model of spins are calculated as gradients of smoothly-varying continuous functions [51].
More generally, however, the spin-lattice interaction should be expanded in terms that are bilinear in the spins S⃗i and S⃗j of atoms i and j, including i = j, and with the atomic coordinates and vector derivatives of J ij , with respect to the coordinates, occurring in increasing order (See Ref. [54]). In such a scheme, the generalised Heisenberg exchange term in Eq. (3.16) corresponds to the zeroth-order term in the bilinear spin expansion. The next-lowest order term, the first-order term of the bilinear expansion, which is not included in Eq. (3.16), contains a sum over all atoms of the dot product of a rank 3 tensor, the gradient of J ij with respect to the position of atom k, and the position vector of atom k [54].
Interestingly, the (anisotropic) off-diagonal elements of the rank 3 tensor in the aforementioned first-order term of the bilinear expansion provides a (non-relativistic) means for the lattice and spins to exchange angular momentum and equilibrate to the same temperature in a simulation in the microcanonical ensemble [54–56]. In the next chapter, it will be shown that the transfer of angular momentum between spin and lattice degrees of freedom, the Einstein-de Haas effect [51,163], can be achieved in SLD simulations of ferromagnetic nanocontacts via the addition of a generalised uniaxial magnetocrystalline anisotropy correction [55,56].
The symmetry properties of the tensor in the first-order term, as well as those of higher rank in the higher order terms, of the bilinear expansion, depend crucially on the geometry and symmetry of the systems under consideration. Hence, in Ref. [54], the elements of the rank 3 tensor in the first-order term were only worked out for a dimer, trimer and tetramer of iron atoms, in addition to bulk iron. Repeating such an analysis for an iron nanocontact of arbitrary shape, which evolves dynamically as it stretched or compressed, is well outside the scope of this thesis. Besides, Wang et al. showed that although the model of J ij in Eq. (3.17) is too simple [53], J ij (rij ) could be represented by a superposition of (isotropic)
Bethe-Slater curves.
Finally, returning to Eq. (3.16), note that a term J ij (rij )(1) has been subtracted from the spin-dependent generalised Heisenberg exchange term to ensure that the energy and forces are properly defined in the collinear ferromagnetic phase that exists at 0 K [48]. This term represents the ground-state energy of the spin degrees of freedom at 0 K. Additionally, it permits the use of any EAM interatomic potential, not restricted to the Derlet and Dudarev potential.
In this regard, because the EAM potential developed by Malerba et al. in Ref.
[43] reproduces the energies of exposed (001), (110) and (111) surfaces of Fe better than any of the other potentials, including the Derlet-Dudarev potential, it will be used in the production runs carried out in this work.
The equations of motion integrated during simulations in SPILADY are, written succinctly as, It is noted that the sign in the Hamilton’s equation, Eq. (3.18), for the force contribution ⃗fi =−∇r⃗i Η from the generalised Heisenberg exchange term in Eq. (3.16) is wrong in many of Ma et al.’s works (e.g., Refs. [48,134]). The equations are, however, correctly expressed in e.g., Refs. [49,50,164]. Notwithstanding, the sign of this force contribution has been verified and found to be correct in the SPILADY code itself.
Technical details of the integration procedure and temperature control in SLD simulations is explained at length in Refs. [48,134,165,166]. A more sophisticated integration scheme, the second-order Suzuki-Trotter decomposition (STD) of non-commuting operators of spin and lattice coordinates (see Ref. [48,166–168] and references therein), simultaneously conserving the total energy, linear momentum and the spin magnitudes of the atoms for a reasonable trade-off between accuracy and computational efficiency [48,166–168], is used in SPILADY. The default simulation time step in SPILADY is one femtosecond, the same as in LAMMPS.
Finally, the temperature in SPILADY is controlled by a Langevin thermostat (for a very detailed discussion of the implementation thereof, see Ref. [165]). If so desired, two separate thermostats can be used simultaneously to thermalise the lattice and spin degrees of freedom. Further details regarding temperature control will be provided in the next chapter, which deals with adding the effect of spin-orbit coupling to SLD simulations. Thermostatting in LAMMPS is discussed in Chapter 5.
In this work, two different open-source codes are used to model metallic nanocontacts. The first, the Large-scale atomic/molecular massively parallel simulator (LAMMPS) [169], is a well-established and versatile simulator, which has been used extensively, in combination with EAM potentials, throughout our previous work [23,24,76].
Spin-lattice dynamics, on the other hand, is still in its infancy. So far, it has only been applied to a few, selected, systems [48,49,55,134]. In the present work, SLD is applied for the first time to ferromagnetic nanocontacts. For this reason, a significant portion of the present work is devoted to the extension of the SLD model and validation thereof, as will be explained in detail in the next chapter.

CONTENTS
Declaration
Dedication
Acknowledgments
Abstract
1. Introduction
1.1. Background
1.2. Thesis statement and research objectives
1.3. Approximations and limitations
1.4. Definitions of technical terms and abbreviations
1.5. Significance
2. Literature review
2.1. A brief account of metallic nanocontacts and scalar-relativistic effects
2.2. An account of domain-wall magnetoresistance in Ni and Fe nanocontacts
3. Overview of standard methods
3.1. Classical molecular dynamics
3.2. Spin-lattice dynamics
3.3. Density functional theory
3.4. Non-equilibrium Green’s Function DFT quantum transport
4. Extensions of standard methods and preliminary computational results
4.1. Emulating the experiments: cyclic loading
4.2. Spin-lattice dynamics with spin-orbit coupling
4.3. Vector-relativistic NEGF quantum transport
5. Relativistic effects in non-magnetic metal nanocontacts: Au, Cu and Ag
5.1. Classical molecular dynamics of Au, Ag and Cu nanocontacts
5.2. The role of scalar-relativistic effects in Au
6. The atomic configurations of Ni and Fe before rupture
6.1. Why theoretical conductance values for BCC Fe nanocontacts are so low
6.2. Comparison with an FCC ferromagnetic metal: nickel
7. The role of ferromagnetism and uniaxial magnetic anisotropy in Fe and Ni nancontacts
7.1. Magnetic domain walls in (001)- and (111)-oriented Ni nanocontacts
7.2. Magnetic domain walls in (001)-oriented Fe nanocontacts
8. Concluding remarks
Bibliography
GET THE COMPLETE PROJECT