Signatures of electron-magnon interaction in charge and spin currents through magnetic tunnel junctions: A nonequilibrium many-body perturbation theory approach

Signatures of electron-magnon interaction in charge and spin currents through magnetic tunnel junctions: A nonequilibrium many-body perturbation theory approach

Farzad Mahfouzi mahfouzi@udel.edu Department of Physics and Astronomy, University of Delaware, Newark, DE 19716-2570, USA    Branislav K. Nikolić bnikolic@udel.edu Department of Physics and Astronomy, University of Delaware, Newark, DE 19716-2570, USA
Abstract

We develop a numerically exact scheme for resumming certain classes of Feynman diagrams in the self-consistent perturbative expansion for the electron and magnon self-energies in the nonequilibrium Green function formalism applied to a coupled electron-magnon (e-m) system driven out of equilibrium by the applied finite bias voltage. Our scheme operates with the electronic and magnonic GFs and the corresponding self-energies viewed as matrices in the Keldysh space, rather than conventionally extracting their retarded and lesser components, which greatly simplifies translation of diagrams into compact mathematical expressions and their computational implementation. This is employed to understand the effect of inelastic e-m scattering on charge and spin current vs. bias voltage in F/I/F (F-ferromagnet; I-insulating barrier) magnetic tunnel junctions (MTJs), which are modeled on a quasi-one-dimensional (quasi-1D) tight-binding lattice for the electronic subsystem and quasi-1D Heisenberg model for the magnonic subsystem. For this purpose, we evaluate Fock diagram for the electronic self-energy and the electron-hole polarization bubble diagram for the magnonic self-energy. The respective electronic and magnonic GF lines within these diagrams are the fully interacting ones, thereby requiring to solve the ensuing coupled system of nonlinear integral equations self-consistently. Despite using the quasi-1D model and treating e-m interaction in many-body fashion only within a small active region consisting of few lattice sites around the F/I interface, our analysis captures essential features of the so-called zero-bias anomaly observed [Phys. Rev. B 77, 014440 (2008)] in both MgO- and AlO-based realistic 3D MTJs where the second derivative (i.e., inelastic electron tunneling spectrum) of charge current exhibits sharp peaks of opposite sign on either side . We show that this is closely related to substantially modified magnonic density of states (DOS) after e-m interaction is turned on—the magnonic bandwidth over which DOS is non-zero becomes broadened, thereby making e-m scattering at arbitrary small bias voltage possible, while DOS also acquires peaks (on the top of a continuous background) signifying the formation of quasibound states of magnons dressed by the cloud of electron-hole pair excitations. We also demonstrate that the sum of electronic spin currents in all of the semi-infinite leads attached to the active region quantifies the loss of spin angular momentum carried away from the active region by the magnonic spin current.

pacs:
72.25.-b, 72.10.Di, 72.10.Bg

I Introduction

Magnetic tunnel junctions (MTJ) are layered heterostructures in which an insulating tunnel barrier (I) separates two ferromagnetic layers (F). They have been the subject of vigorous research in both fundamental and applied physics since they exhibit effects like tunneling magneto-resistance (TMR) Tsymbal2003 () and spin-transfer torque (STT), Ralph2008 (); Brataas2012 () as well as quantum size effects in electron transport (even at room temperature) when normal metal (N) layer is inserted. Yuasa2002 () From the fundamental viewpoint, these effects represent examples of nonequilibrium quantum many-body systems with an interplay of fast conduction electrons carrying spin current and slow collective magnetization, while from the viewpoint of applications they play an essential role in developing magnetic sensors, random access memory, novel programmable logic devices, resonant-tunneling spin transistors and nanoscale microwave oscillators with ultrawide operating frequency ranges. Katine2008 ()

The STT is a phenomenon in which a spin current of sufficiently large density injected into F layer either switches its magnetization from one static configuration to another or generates a dynamical situation with steady-state precessing magnetization.Ralph2008 (); Brataas2012 () The origin of STT is the absorption of the itinerant flow of angular momentum components normal to the magnetization direction. For F/I/F MTJs illustrated in Fig. 1, where the reference F layer with fixed magnetization plays the role of an external spin-polarizer and right F layer has free magnetization , it is customary to analyze the parallel and perpendicular components of the STT vector . These two terms act differently on magnetization dynamics— effectively changes the magnetic damping (i.e., an antidamping or additional damping depending on the current polarity), whereas acts like a magnetic field. Ralph2008 ()

A majority of theoretical studies of TMR or STT effects has assumed phase-coherent tunneling of non-interacting quasiparticles. For example, such approaches Butler2001 (); Mathon2001 () have led to a remarkable prediction of very large TMR ratio % at zero bias voltage for clean epitaxial MgO-based MTJs. The TMR ratio is defined by , where is the resistance for parallel orientations of two magnetizations in F/I/F MTJ and is the resistance when they are antiparallel. These predictions have ignited large experimental efforts that have eventually reached TMR ratios of more than 1000% at low temperatures and % at room temperature for well-oriented MgO barriers with stress relaxation. Ikeda2008 () The phase-coherent calculations—such as numerical ones based on the nonequilibrium Green function (NEGF) formalism Stefanucci2013 () combined with simplistic tight-binding Hamiltonians Theodonis2006 (); Tang2009 () or first-principles obtained Hamiltonians; Heiliger2008 () as well as analytical ones Xiao2008a () based on the scattering approach—have been able to capture the dependence of and on the bias voltages V in MgO-based MTJs. Kubota2008 () However, such theories Theodonis2006 (); Tang2009 (); Xiao2008a () including only elastic electron tunneling start to deviate from experimental findings at higher bias voltages, which is particularly pronounced Li2008b (); Park2011a (); Wang2011 () for (playing a significant role during magnetization switching at V).

The inelastic electron-magnon (e-m) or electron-phonon (e-ph) scattering could account for these discrepancies. Li2008b (); Park2011a () In particular, since magnon bandwidth is usually of the order of meV, at high bias voltages multiple magnon scattering events can be excited. Balashov2008 () Also, energy dependence of the magnon density of states (DOS) probed Balashov2008 () at finite bias voltage is intimately linked to the evolution of the magnetization during current-driven switching when going beyond the macrospin approximation. Strachan2008 (); Brataas2006a ()

Figure 1: (Color online) Schematic view of a quasi-1D model of F/I/F MTJ where the left semi-infinite ideal F lead, modeled as spin-split tight-binding lattice of size , is attached via hopping to an active region consisting of lattice sites (we use and or in the calculations below). The right semi-infinite lead is attached to the active region via smaller hopping that simulates the tunnel barrier . The same sites also host localized spins which are coupled to each other via the ferromagnetic coupling in the Heisenberg model. The local coupling between the spin of conduction electrons and localized spin on each site is of strength . The left and right semi-infinite leads are assumed to terminate into macroscopic Fermi liquid reservoirs held at electrochemical potentials and , respectively, whose difference sets the bias voltage . The voltage profile across MTJ is shown on the top. The left F layer is assumed to be attached at infinity to a macroscopic reservoir of magnons held at temperature .

In fact, even at small bias voltage thermally excited magnons affect TMR (e.g., emission or absorption of magnon at F/I interface reduces the effective spin polarization of electrons incoming from F leads in Fig. 1), so that TMR decreases with increasing temperature. Drewello2008 (); Khan2010 () Although thermally induced change of the resistance is different for AlO- and MgO-based MTJs, their inelastic tunneling spectra Reed2008 () (IETS) shows very similar properties. That is, plotting the second derivative of current vs. bias voltage in MTJs reveals zero bias anomaly (ZBA) where peaks (see, e.g., Fig. 2 in Ref. Drewello2008, ) of opposite sign appear at mV and are related to magnons. Also, additional phonon peaks are found Drewello2008 () at mV for the MgO-based MTJs or at mV for AlO-based MTJs.

Theoretical efforts to capture e-m inelastic scattering effects on TMR, ZBA and STT have thus far utilized simplified frameworks Levy2006 (); Manchon2009a () which cannot deal with multiple scattering events, backaction of magnons driven far from equilibrium and energy dependence of magnonic DOS. Such effects can be taken systematically and rigorously into account at arbitrary temperature or bias voltage by using the NEGF formalism coupled with perturbation expansion of electron or magnon self-energies in the presence of their mutual interaction in terms of the respective Feynman diagrams. Stefanucci2013 () In equilibrium problems, like that of magnetic polaron, electronic self-energy has been constructed by considering large set of diagrams involving an arbitrary number of e-m scattering vertices between the emission and absorption vertices. Richmond1970 () However, using the same diagrams within the NEGF framework would violate charge conservation, yielding different charge currents in the left and right lead of a two-terminal device at finite bias voltage.

One of the conserving approximations is the so-called self-consistent Born approximation (SCBA) where one considers Hartree and Fock diagrams for the electronic self-energy which corresponds to perturbation in the order where is the strength of e-m interaction. Evaluation of these diagrams for a systems defined on the lattice hosting orbitals in real space is computationally very demanding due to the fact that GF lines in the diagrams of nonequilibrium many-body perturbation theory (MBPT) are fully interacting (or dressed), so that self-energy matrix becomes a functional of the GF matrix. This generates a coupled system of nonlinear integral equations which has to be solved by performing multiple integrations over each matrix element until the self-consistency is achieved. Such route has been undertaken in only a handful of studies where further simplifications (such as using dispersionless magnons, ) were utilized. Reininghaus2006 ()

Here we discuss in Sec. III how to construct the electronic self-energy and GF within SCBA, together with the magnonic self-energy within the electron-hole (e-h) polarization bubble approximation which takes into account influence of electrons on magnons while inserting the dressed magnonic GF into SCBA diagrams. Thus, consideration of such diagrams is akin to the self-consistent GW treatment of the one-particle electronic self-energy due to electron-electron interaction out of equilibrium. Thygesen2008 () Our approach treats these quantities as matrices Kita2010 () in the Keldysh space, rather than following the commonly used route based on Langreth rules to manipulate expressions involving products of their submatrices. Stefanucci2013 () This formalism is then applied to a many-body Hamiltonian, introduced in Sec. II, of an interacting e-m system defined on the real-space lattice describing MTJ that is brought out of equilibrium by the applied finite bias voltage. Despite using a quasi-one-dimensional (quasi-1D) model for MTJ illustrated in Fig. 1, where e-m interaction is treated diagrammatically only within few lattice sites (denoted as “active region” in Fig. 1) of the left F layers, charge current versus bias voltage and its second derivative obtained in Sec. IV capture essential features of ZBA observed in realistic junctions. The sum of spin currents carried by electrons in the non-interacting F leads attached to this active region is non-zero which, therefore, allows us to quantify in Sec. IV the amount of lost angular momentum of electronic subsystem that is carried away by magnonic spin current. We conclude in Sec. V.

We also provide two Appendices. Appendix A proves that charge current is conserved when electronic GF is computed by including the Fock diagram only, on the proviso that its electronic GF line is dressed (thereby requiring self-consistency) while its magnonic GF line can be bare or dressed. Appendix B shows numerical implementation of the Hilbert transform, utilized to obtain results in Sec. IV, for a function computed on the mesh of adaptively selected energy points which can be arbitrarily spaced form each other.

Ii Hamiltonian for coupled electron-magnon system within MTJ

To make the discussion transparent, we focus on the particular example of e-m interacting many-body system out of equilibrium which emerges within the quasi-1D model of a two-terminal F/I/F MTJ depicted in Fig. 1. This can be described by the following Hamiltonian

(1)

The first term in Eq. (II) accounts for the on-site potential due to the voltage profile shown in Fig. 1, as well as for the coupling of itinerant electrons to collective magnetization described by the material-dependent exchange potential eV. We use the standard notation for the vector of the Pauli matrices, where are their matrix elements, and is the unit matrix. Here or depending on whether the magnetization of the left or right F layer is parallel or antiparallel to the -axis, respectively. The second term in Eq. (II) describes hopping of electrons between single -orbitals located on the tight-binding lattice sites, where () creates (annihilates) electron on site in spin state and is the nearest-neighbor hopping parameter. We set eV for all pairs of lattice sites, except for the last row of sites of the left F layer and the first row of sites of the right F layer where eV (in the case of active region in Fig. 1) or eV (in the case of active region in Fig. 1) simulates the presence of the tunnel barrier .

The third term is the Heisenberg model Rastelli2013 () describing interaction between spin operators and localized on the nearest-neighbor sites of the same tight-binding lattice, where ferromagnetic coupling is set as eV, except for the I region where . The fourth term with eV is introduced to select energetically favorable direction (i.e., an easy-axis) for the spontaneous magnetization in the ferromagnetic layers along the -axis. Finally, the fifth term describes interaction of the spin operator of conduction electrons with the localized spin operators , where the coupling constant is set as eV.

The active region of MTJ in Fig. 1, within which NEGFs and self-energies due to e-m interaction are computed, consists of sites enclosed in Fig. 1. The rest of the tight-binding sites belong to the left and right semi-infinite leads (taken into account through lead self-energies discussed in Sec. III). The leads are assumed to terminate at infinity into macroscopic Fermi liquid reservoirs held at electrochemical potentials and (the Fermi energy is chosen as eV for active region and eV for active region), whose difference sets the bias voltage . Concurrently, the left F layer is assumed to be attached at infinity to a macroscopic reservoir of magnons held at temperature .

Using the approximate version of the Holstein-Primakoff transformations Rastelli2013 ()

(2a)
(2b)
(2c)

we can replace the spin operators by bosonic operators. The approximation in Eq. (2) is a valid when the occupation number of bosonic states at temperature is low, , where we select . This is equivalent to saying that the left or right F layer is near its ferromagnetic ground state where on all lattice sites within both the left and right F layers for the parallel (P) configuration of magnetizations in MTJ, or within the left F layer and within the right F layer for antiparallel (AP) configuration of magnetizations in MTJ.

This replacement allows us to rewrite the Hamiltonian in Eq. (II) as , where all three terms are now given in the second quantization

(3a)
(3b)
(3c)

Here is the number of nearest neighbor sites. The many-body interaction is encoded by in Eq. (3c), which is assumed to be non-zero only in the active device region in Fig. 1. Its first term has a clear physical interpretation—the spin of a conduction electron is flipped when magnon is absorbed or emitted. Since its second term (in the lowest order, its role is to renormalize the effective Zeeman splitting for electrons) is much smaller than the first term due to assumed , we ignore it in subsequent discussion.

In the rest of the system depicted in Fig. 1, electrons and magnons are assumed to behave as non-interacting quasiparticles. For example, in the case of active region, the left and right semi-infinite F leads are 1D chains whose electrons in equilibrium () are described by in Eq. (3a) generating dispersion which is spin split both by the mean-field treatment of interaction term and . The left semi-infinite lead Meier2003 () for the magnonic subsystem is described by in Eq. (3b), where non-interacting magnons have dispersion .

Iii Nonequilibrium diagrammatics for electron-magnon interacting system

iii.1 Compact analytical expressions in the Keldysh space

The NEGF formalism Stefanucci2013 () operates with two central one-particle quantities—the retarded GF ( for fermions or for bosons), describing the density of available quantum states; and the lesser GF ( for fermions or for bosons), describing how quasiparticles occupy those states. One can also use two additional GFs, advanced [ for fermions or for bosons] and greater ( for fermions or for bosons), describing the properties of the corresponding empty states. These four GFs, which generally depend on two time arguments , are connected by the fundamental relation

(4)

for electronic () or bosonic () GFs.

Besides having a clear physical meaning, these four GFs make it possible to obtain nonequilibrium expectation values of any one-particle observable, such as charge and spin currents that are the focus of our study. However, these four GFs do not have perturbation expansion akin to zero-temperature GFs on the Feynman contour (the real-time axis from to ) or finite-temperature GFs on the Matsubara contour (a segment on the imaginary-time axis from to , where ). Instead, perturbation expansion is formulated for the contour-ordered Stefanucci2013 () GF whose two time arguments are located on the Keldysh-Schwinger contour consisting of two counter-propagating copies of the real-time axis—the forward branch extending from to and the backward branch extending from to . Equivalently, one can introduce matrix GF in the so-called Keldysh space Kita2010 () which depends on the two time arguments located on the single real-time axis extending from to . Such Keldysh-space matrix GF for fermions is defined by

(5)

and for bosons it is defined by

(6)

Here and ; is the Heaviside step function ( for and for ); and is the nonequilibrium expectation value Stefanucci2013 () with being the contour ordering operator and the initial density matrix of the system is usually taken at for the steady-state formulations within the NEGF. The Keldysh-space matrices (such as or ) satisfy

(7)

where we use to denote the Pauli matrices acting in the Keldysh space.

In stationary problems GFs depend only on the time difference , so that they can be Fourier transformed to energy or frequency . Using Hamiltonian in Eq. (3) and performing such Fourier transform leads to the following Keldysh-space Dyson equation for electrons

(8)

or magnons

(9)

This approach allows for compact notation by avoiding the widespread route Stefanucci2013 () where one starts from the Dyson equation for the contour-ordered GF containing convolution integrals on the two-branch Keldysh-Schwinger contour, and then applies the so-called Langreth rules Stefanucci2013 () to find the lengthy expressions involving the lesser and retarded GFs with two time arguments located on the single real-time axis (or their Fourier transforms). The Dyson equation for the Keldysh-space matrix GFs in energy, like Eqs. (8) and  (9), is rarely used in the literature due to redundancy expressed by Eq. (4). For example, such equation can be found in the NEGF-based calculations of the full counting statistics Gogolin2006 (); Urban2010 (); Novotny2011 () where the presence of the counting field results in the nonunitary evolution on the Keldysh-Schwinger contour, thereby requiring to work with all four submatrices in Eqs. (5) or  (6) because the relation like Eq. (4) is not valid anymore. Nevertheless, even when Eq. (4) holds, it can be advantageous to work directly with the matrices in the Keldysh space—this greatly simplifies writing the analytical expressions that the Feynman diagrams of nonequilibrium MBPT represent and, moreover, it makes possible to derive expressions for the perturbation expansion of more complicated quantities like the current noise obtained from the two-particle nonequilibrium correlation function. Mahfouzi2013b ()

Figure 2: (Color online) Diagrammatic representation of the Dyson equation in the Keldysh space: (a) the electron case, , in Eq. (8); and (b) the magnon case, , in Eq. (9). The perturbation expansion for the electronic self-energy in (a) retains Hartree and Fock diagrams, while the expansion in (b) for the magnonic self-energy retains e-h polarization bubble diagram. The single straight line denotes the non-interacting electronic GF, [which includes the self-energies due to the leads in Eq. (10)]; the single wavy line denotes the non-interacting magnonic GF, [which includes the self-energy due to the left lead in Eq. (11)]; double straight line denotes the interacting electronic GF, ; and double wavy line denotes the interacting magnonic GF, . The solid circles denote vertices that are integrated out. The electron spin is flipped at each vertex, which is illustrated by spin- (before the vertex) being flipped into spin- (after the vertex), while a magnon is being created. The same process applies to flipping of spin- to spin-, where the direction of magnon propagation (indicated by arrow on the wavy lines) is reversed. Note that the Hartree diagram in panel (a) contains a single (rather than double) wavy line in order to avoid double counting.

The self-energies due to many-body interaction, for electron and for magnon, can be simply added to the self-energies introduced by the attached semi-infinite leads, for electron and for magnons, respectively. This is due to the fact that e-m interaction is assumed to be localized within the active region in Fig. 1, so that the leads do not involve many-body interactions. Thus, the self-energies of the leads for the junction in Fig. 1

(10)
(11)

are single-particle quantities which can always be computed in an exact fashion, either analytically Datta1995 () for simple models like Eqs. (3a) and  (3b) or numerically for more complicated lead Hamiltonians. Rungger2008 () The effect of the bias voltage is introduced by a rigid shift in energy, . The Fermi function of the macroscopic reservoir to which lead is attached is denoted by . The Bose-Einstein distribution function of the macroscopic reservoir of magnons to which the left F layer is attached is denoted by .

The one-particle self-energies due to e-m interaction, and magnon , are formally obtained by summing all irreducible diagrams, i.e., those diagrams that cannot be taken apart by cutting a single line. The self-energies are actually functionals of the respective electronic or magnonic GF, so that they have to be approximated in practical calculations. Diagrammatic techniques provide a natural scheme for generating approximate self-energies, as well as for systematically improving these approximations. While there are no general prescriptions on how to select the relevant diagrams, this process can be guided by physical intuition. In addition, unlike equilibrium Richmond1970 () MBPT, the diagrams selected in nonequilibrium Stefanucci2013 () MBPT must generate GFs that yield expectation value for charge current which is conserved. For example, the final current in the left and right leads of the device in Fig. 1 must be the same in any chosen approximation for the self-energies.

Here we select one of the conserving approximations, Stefanucci2013 () where the Feynman diagrams retained for the electron or magnon self-energies are displayed in Figs. 2(a) and  2(b), respectively. The diagrams in Figs. 2(a) are equivalent to the so-called self-consistent Born approximation (SCBA) for the electronic self-energy considered in problems like e-ph interacting systems. Mitra2004 (); Viljas2005 (); Koch2006 (); Frederiksen2007 (); Galperin2007 (); Lu2007 (); Asai2008 (); Lee2009a (); Dash2011 () However, in the case of e-m interaction, one has to introduce additional bookkeeping in such diagrams to account for flipping of electron spin together with magnon emission or absorption, as illustrated in Fig. 2. The Keldysh-space expression for the electronic self-energy read from the Fock diagram in Fig. 2(a) is given by

(12a)
(12b)

Here and are matrix indices which include the Keldysh space and real (i.e., orbital) space, so that selects a submatrix of .

Finding the proper expression for the Hartree diagram in Fig. 2(a) in the Keldysh space requires extra care Oguri2013 () because and for the inner electronic GF along the loop are equal, so that in Eq. (5) gives 1 or 0 depending on whether or , respectively. Therefore, in terms of which we obtain the following expressions

(13a)
(13b)

Note that the off-diagonal (i.e., lesser and greater) components of vanish, and the remaining retarded component on the diagonal is energy independent.

Although in SCBA, we retain only the Fock term in the actual calculations below. We note that in SCBA for e-ph interacting systems, is often neglected Frederiksen2007 (); Lee2009a () due to being small and, therefore, having little effect on the final current (this becomes unwarranted for larger e-ph interaction strengths where SCBA breaks down Lee2009a ()). For the e-m interacting systems, the situation is much more complex because direct evaluation of Eq. (13) leads to numerical instabilities. This stems from the fact that our MTJ is invariant with respect to the rotation around the -axis (see Fig. 1), so that spin-flip rate which appears in Eq. (13) can acquire arbitrary phase thereby requiring to consider full double time dependence of . We relegate this to future studies, while here we retain which is termed Lee2009a () Fock-only SCBA (F-SCBA).

In the case of e-ph many-body systems driven far from equilibrium, phonon heating due to propagating electrons has been considered either phenomenologically using a rate equation for the phonon occupation, Koch2006 (); Frederiksen2007 (); Siddiqui2007 () or microscopically by using phonon GF with interacting self-energy truncated to the e-h polarization bubble diagram. Mitra2004 (); Viljas2005 (); Galperin2007 (); Lu2007 (); Asai2008 (); Urban2010 (); Novotny2011 () It is worth mentioning that the two approaches yield virtually identical results for time-averaged current in the limit of weak e-ph coupling, but they start differing significantly in the case of the current noise due to the feedback of the phonon dynamics on the statistics of the transmitted electrons which cannot be captured by the phenomenological rate equation approach. Novotny2011 () Since magnon bandwidth ( meV) is relatively small, Balashov2008 () they can be easily driven into far from equilibrium state by charge current at finite bias voltage. For the purpose of describing such state, we retain in Fig. 2(b) the e-h polarization bubble diagram for the magnonic self-energy whose analytical expression is given by

(14)

Thus, the dressed magnonic GF which includes this self-energy through Eq. (9) will be inserted into the electronic self-energy diagrams in Fig. 2(a), thereby generating an infinite resummation of diagrams until the mutual self-consistency is achieved. Note that this is analogous to the self-consistent GW treatment of the one-particle electronic self-energy due to electron-electron interaction out of equilibrium. Thygesen2008 ()

iii.2 Numerical implementation in Keldysh space

Equations (8),  (9),  (12) and  (14) form a system of coupled nonlinear integral equations that has to be solved iteratively until the convergence criterion is met. We use expectation value of charge current (see Sec. IV) to define one such criterion, . Here is charge current in lead at the beginning of an iteration, denotes charge current at the end of the same iteration and we select .

Figure 3: (Color online) Elastic and inelastic contributions to charge current in the left and right lead of MTJ in Fig. 1 with active region at different iteration number within the self-consistent loop for solving coupled Eqs. (8),  (9),  (12) and  (14). The orientation of the magnetizations of two F layers in Fig. 1 is parallel in (a) and antiparallel in (b). The bias voltage is set as mV and temperature is K.

The Keldysh-space electronic GF and self-energies in these coupled equations are matrices of the size (if orbitals are used per site then ), while the magnonic GF and self-energies are matrices of the size . The most time-consuming part of solving the coupled equations is the integration in Eq. (12) for , which can be viewed as the convolution

(15)

where is the elementwise product of matrices. The fact that matrix elements of the retarded and advanced components of GFs in Eqs. (8) and  (9) are nonzero in the whole range of integration would make the numerical computation of this convolution prohibitively expensive. However, this obstacle can be removed by using the fact that retarded and advanced GFs are analytic functions in the upper and lower half of the complex plane, respectively. Thus, the real () and imaginary () parts of their matrix elements must obey the following relation

(16)

Here is the Hilbert transform, whose implementation in our scheme is discussed in more details in Appendix B, and stands for the Cauchy principal value. This makes it possible to decompose Keldysh-space matrices as follows

(17)

where

(18a)
(18b)

are labeled as “symmetric” and “asymmetric” part. We note that all the relevant information is already contained in , so that one can find from it by using

(19)

This idea allows us to restrict the range of integration in the convolution in Eq. (15) to the energy bandwidth of electrons and magnons. Once the decomposition is done for both matrices and , the convolution in Eq. (15) can be evaluated using

(20)

with

(21a)
(21b)

Here we used the following properties of the Hilbert transform and the convolution operator

(22)
(23)

Note that one has to actually calculate only , after which the asymmetric part is obtained from Eq. (19).

Figure 4: (Color online) Total spin current dissipated (at K) inside active region of MTJ in Fig. 1 as a function of the bias voltage for parallel and antiparallel orientation of the magnetizations of two F layers. The spin current is obtained from Eq. (28) using electronic GF computed by solving coupled Eqs. (8),  (9) and  (12).

Iv Application to charge and spin currents in MTJ driven by finite bias voltage

The charge current in lead can be viewed as the sum of two spin-resolved charge currents, . For an interacting active region attached to two non-interacting semi-infinite leads, lead currents can be obtained directly from and

(24)

where we employ the following notation

(25)

The second line in Eq. (IV) is the well-known Meir-Wingreen formula. Meir1992 () Similarly, the spin current (in the same units as for the charge current) in lead is obtained from

(26)
Figure 5: (Color online) The electronic density of states within the active region of MTJ model in Fig. 1 of size: (a) ; and (b) . The magnonic density of states within the same active regions is shown in panels (c) and (d), respectively. These quantities are computed at finite bias voltage mV and at temperature K, in the absence ( for dashed line) or the presence ( for solid line) of e-m interaction in the Hamiltonian in Eq. (3c). The respective DOS is obtained from the retarded electronic GF, using , or the retarded magnonic GF, using , after solving coupled Eqs. (8),  (9),  (12) and  (14) which take into account influence of electrons on magnons. The arrows in panels (a) and (b) point at the kinks (located at the Fermi energies of MTJ model with two different sizes of the active region, respectively) in the interacting electronic DOS (solid line) due to e-m coupling.
Figure 6: (Color online) Elastic and inelastic charge currents in Eq. (27) and their sum, as well as the corresponding first and second derivatives, versus the bias voltage in the model of MTJ in Fig. 1 with active region for parallel and antiparallel orientation of its magnetizations. These charge currents are obtained from the electronic GF computed in F-SCBA which includes (in self-consistent fashion) the non-interacting magnonic GF. The temperature is set as K.

The charge current in Eq. (IV) can be conveniently separated Asai2008 (); Ness2010 () into two terms,

(27a)
(27b)

We label the first term as “elastic” current since it has the form of the Landauer-like formula for elastic transport of non-interacting quasiparticles whose effective transmission function is expressed Caroli1971 () in terms of NEGF quantities, . The second term appears as the nonequilibrium corrections due to many-body interactions, which we label as “inelastic” current. Plotting separately elastic and inelastic current components makes it possible to provide additional insights when interpreting our results in Sec. IV.

Figure 7: (Color online) Elastic and inelastic charge currents in Eq. (27) and their sum, as well as the corresponding first and second derivatives, versus the bias voltage in the model of MTJ in Fig. 1 with active region for parallel and antiparallel orientation of its magnetizations. These charge currents are obtained from the electronic GF computed in F-SCBA, which includes (in self-consistent fashion) the interacting magnonic GF with the e-h polarization bubble diagram in Fig. 2(b). The temperature is set as K.

Note that apparent connection of Eq. (27a) to the Landauer formula should not be pushed too far since the effective transmission in already contains part of e-m interaction. That is, the standard Landauer formula Caroli1971 () for single-particle elastic scattering uses the retarded and advanced GFs which include the self-energies due to the semi-infinite leads only. On the other hand, and in include additional self-energy due to e-m interaction which renormalizes the non-interacting reference system, and for strong enough interaction can go even beyond the quasiparticle description of the many-body interacting quantum system. Even when quasiparticles are well-defined, the presence of self-energy that is functional of the retarded GF itself means that includes dephasing effects due to many-body interaction Chen2012 () and is, therefore, different from phase-coherent tunneling current that would be obtained from the the standard Landauer formula. Caroli1971 ()

Here we illustrate in Fig. 3 that is conserved at each iteration, while the conservation of component requires to reach the self-consistency in the computation of electronic GF and self-energy in F-SCBA, as discussed in Eq. (A). Note that the magnonic GF and self-energy used to obtain Fig. 3 also include e-h polarization bubble diagram from Fig. 2(b).

The spin current can analogously be separated into the elastic and inelastic contributions,

(28a)
(28b)

We find that vanishes at all bias voltages, so that the total spin current is governed by the inelastic component only which is plotted in Fig. 4. Thus, this quantity measures the loss of angular momentum of electrons within the interacting active region of MTJ in Fig. 1, which is then carried away by magnonic spin current (through the left semi-infinite lead toward the left magnonic macroscopic reservoir). Although spin current carried by electrons or magnons individually is not conserved, the total angular moment in this process is conserved. This finding further justifies the separation of currents into elastic and inelastic contributions since does not participate in the loss of angular momentum.

Figure 8: (Color online) Elastic and inelastic charge currents in Eq. (27) and their sum, as well as the corresponding first and second derivatives, versus the bias voltage in the model of MTJ in Fig. 1 with active region for parallel and antiparallel orientation of its magnetizations. These charge currents are obtained from the electronic GF computed in F-SCBA, which includes (in self-consistent fashion) the non-interacting magnonic GF. The temperature is set as K.

Due to the fact that the e-m interaction strength is comparable to the magnonic bandwidth, the single particle and many-body properties of magnons within the active region in Fig. 1 are governed largely by the collective quasiparticles rather than the bare (non-interacting) magnons we started from. This is demonstrated by plotting the magnon density of states (DOS) in Fig. 5 within the active region versus energy. The DOS is obtained from with e-m interactions turned off () or turned on (). In Fig. 5(c), we can clearly distinguish three peaks corresponding to the quasibound states suggesting the formation of long-lived quasiparticles. They can be interpreted as a magnon dressed by the cloud of electron-hole pair excitations out of equilibrium. Importantly for ZBA discussed below, the DOS of interacting magnons extends all the way to zero energy, thereby enabling e-m scattering even at vanishingly small bias voltage. On the other hand, the electronic DOS obtained from and plotted in Figs. 5(a) and  5(b) is only slightly perturbed when e-m interaction is turned on due to much larger electronic bandwidth.

Figure 9: (Color online) The TMR vs. bias voltage in the model of MTJ in Fig. 1 with active region . The inset shows TMR as a function of temperature in the linear-response limit . These results are obtained from electronic GF computed by solving coupled Eqs. (8),  (9) and  (12).

Figure 6 plots the elastic, inelastic and total charge currents, together with their first derivative (i.e., differential conductance) and second derivative (i.e., IETS Reed2008 ()), as a function of the applied bias voltage . The currents in Fig. 6 are computed using F-SCBA for the electronic GF and self-energy of active region in Fig. 1, while the magnonic GF is used as the non-interacting one by setting in Eq. (9). The inelastic current in Fig. 6(b) is zero until the threshold bias voltage is reached ( mV according to dashed line in Fig. 5) at which magnons can be excited. Above the threshold voltage, inelastic current displays Ohmic behavior. This is simply due to the fact that the rate of energy (and angular momentum) loss is proportional to the rate of electrons being injected into the active region. Although elastic current in Fig. 6(a) shows apparent Ohmic behavior for all bias voltages, the corresponding differential conductance in Fig. 6(d) deviates strongly from the straight line within the energy rage where magnons can be excited. This can be explained by the fact that the effective electronic DOS inside the active region can be changed through e-m scattering. The elastic differential conductance in Fig. 6(d) decreases once the magnons are excited, but this is compensated by the increase of inelastic differential conductance in Fig. 6(e) such that the total differential conductance in Fig. 6(f) has less pronounced features. The IETS plotted in Figs. 6(g) and 6(h) shows a clear signature Reed2008 () of the non-interacting magnonic DOS from Fig. 5 with two peaks emerging slightly away from . Nevertheless, when these two contributions are summed up in Fig. 6(i), IETS for total current shows more than just two peaks.

In order to see the effect of DOS of interacting magnons (solid line in Fig. 5), or possible magnon heating due to tunneling electrons, Fig. 7 presents the same information as in Fig. 6 but recomputed by including e-h polarization bubble diagram in Fig. 2(b) for the magnonic self-energy . Since DOS of interacting magnons in Fig. 5 is sufficiently broadened to reach low frequencies, the inelastic current in Fig. 7(b) is now non-zero even for very small bias voltage . The presence of magnons dressed by the cloud of e-h pair excitations in this calculation forces elastic conductance in Fig. 7(d) to increase around , or inelastic conductance in Fig. 7(e) to decrease, which is opposite to the behavior of the same quantities in the case of non-interacting magnons analyzed in Fig. 6. The two peaks of opposite sign around in partial IETS plotted in Figs. 7(g) and  7(h) look very similar to ZBA peaks observed experimentally Drewello2008 () in realistic MTJs.

Due to opposite effect of e-m interaction on the two partial IETS, their sum in both Figs. 6(i) and  7(i) looses the simple two peak structure around observed experimentally. Drewello2008 () To investigate whether this complexity in the total IETS could be an artifact of 1D nature of MTJ model (with active region attached to 1D leads) considered in Figs. 6 and  7, we recompute the same quantities for the active region in Fig. 8. For this case, both the partial IETS in Figs. 8(g) and  8(h) and the total IETS in Fig. 8(i) exhibit simple two peak structure. However, the two peaks appear slightly away from the zero bias voltage because we do not include (due to substantial computational expense) e-h polarization bubble diagram from Fig. 2(b) which is needed to introduce non-zero magnonic DOS at low energies in Fig. 5(d) enabling e-m scattering at .

The usage of active region introduces non-negligible ratio which can be extracted from Fig. 8(c) as for . Its detailed dependence on plotted in Fig. 9 shows how ZBA vanishes for temperatures K. Also, the TMR ratio (at ) vs. temperature shown in the inset of Fig. 9 agrees with experimentally observed Drewello2008 (); Khan2010 () TMR decrease with increasing temperature.

Considering fully 3D model of MTJs, where additional -point sampling is required for the transverse direction, would require carefully crafted approximations to evade prohibitively expensive five-dimensional integrals in the systems of coupled nonlinear integral equations for the electronic and magnonic GFs. Also, we note that in experiments Drewello2008 () has a dip at , and its absolute value increases with increasing due to opening of new conducting channels by inelastic e-m scattering. This is not seen in Figs. 6(f),  7(f) and  8(f) due to small number of spin-resolved conducting channels (up to six for active region) present in our model of MTJ.

V Concluding remarks

The e-ph interaction in nanostructures driven out of equilibrium, as the example of nonequilibrium electron-boson quantum-many body system, has been amply studied Mitra2004 (); Viljas2005 (); Koch2006 (); Frederiksen2007 (); Galperin2007 (); Lee2009a (); Dash2011 (); Urban2010 (); Novotny2011 () over the past decade using NEGF formalism. This approach, which makes it possible to rigorously model microscopic details of inelastic scattering processes, has been typically implemented using the SCBA diagrams for the electronic self-energy, and sometimes also including e-h polarization bubble diagram for the phonon self-energy in the nonequilibrium MBPT. On the other hand, the same level of description of e-m scattering has received far less attention, Reininghaus2006 () despite its great relevance for a plethora of problems in spintronics. Levy2006 (); Manchon2009a ()

In this study, we have shown how to obtain analytical expressions for SCBA and e-h polarization bubble diagrams describing e-m scattering. This is achieved in a particularly compact form by using matrix GFs in the Keldysh space (which are functions of energy for electrons or frequency for magnons in steady-state nonequilibrium), thereby simplifying tracking of electron spin flips and direction of magnon propagation required to conserve angular momentum at each vertex of the Feynman diagrams. The self-consistent solution of the corresponding system of coupled nonlinear integral equations, which is equivalent to infinite resummation of certain classes of diagrams (akin to the self-consistent GW treatment of the one-particle electronic self-energy due to electron-electron interaction out of equilibrium Thygesen2008 ()), is obtained via several intertwined numerical algorithms that reduce the computational complexity of this task.

Using this framework, we have computed charge and spin currents at finite bias voltage in quasi-1D models of F/I/F MTJ illustrated in Fig. 1. Our key results are summarized as follows: (i) while elastic component of the sum of spin currents in all attached leads is zero at all bias voltages, the inelastic one is non-zero thereby measuring the loss of spin angular momentum carried by magnons away from the active region (see Fig. 4); (ii) turning on the e-m interaction strongly modifies magnonic DOS, which acquires larger bandwidth while exhibiting peaks due to quasibound states of magnons dressed by the cloud of electron-hole pair excitations [see Figs. 5(c) and  5(d)]; (iii) using F-SCBA for the electronic self-energy in Fig. 2(a), coupled with e-h polarization bubble diagram for the magnonic self-energy in Fig. 2(b), generates two peak structure around zero bias voltage in the second derivative (i.e., IETS Reed2008 ()) of both elastic and inelastic charge currents (see Fig. 7).

We emphasize that e-h polarization bubble diagram in Fig. 2(b) is responsible for the substantial change of magnonic DOS (encoded by the retarded magnonic GF) in equilibrium, as well as for magnon heating (encoded by the lesser magnonic GF) in nonequilibrium due to tunneling electrons where heated magnons can also exert backaction Novotny2011 () onto electrons. Since ZBA occurs at very small bias voltages, the former effect is more important because the broadened magnonic DOS in Figs. 5(c) and  5(d) extends to low energies (in contrast to DOS of non-interacting magnons), thereby making possible inelastic e-m scattering even at zero bias. It is worth mentioning that the effect of electron-boson interaction on bosonic DOS has been rarely discussed in prior NEGF studies Reininghaus2006 (); Mitra2004 (); Viljas2005 (); Koch2006 (); Frederiksen2007 (); Galperin2007 (); Lee2009a (); Dash2011 (); Urban2010 (); Novotny2011 () of coupled electron-boson systems, either due to simplicity of bosonic spectrum assumed or because of not inserting dressed magnonic GF lines (which include e-h polarization bubble diagram) into SCBA for electronic GF.

While partial IETS in Figs. 7(g) and  7(h) obtained from elastic or inelastic components of charge current, respectively, are quite similar to experimentally observed Drewello2008 () ZBA, their sum in Fig. 7(i) is more complicated due to usage of strictly 1D model with active region in those Figures. Switching to active region attached to quasi-1D leads (supporting more than two spin-resolved conducting channels) makes total IETS in Fig. 8(i) exhibiting only two peaks, albeit shifted slightly away from due to exclusion of e-h polarization bubble diagram in this calculation in order to reduce computational expense.

We believe that extension of our approach to 3D MTJs [by adding computationally expensive -point sampling in the and directions while keeping real space Hamiltonian from Eq. (II) in the direction] would be able to describe not just ZBA in realistic junctions, but also TMR and STT effects as a function of temperature and bias voltage, thereby opening a path to understand how to optimize these effects for applications in spintronics by tailoring magnon spectrum.

Acknowledgements.
We thank H. Mera for immensely valuable discussions and comments on several drafts. This work was supported in part by NSF under Grant No. ECCS 1202069.

Appendix A Conservation of charge current in F-SCBA

It is instructive to check if charge current is conserved, , after electronic GF in F-SCBA is inserted into Eq. (IV)

(29)

We use to write the second line in Eq. (A). To show that the third line is identically zero, we use the fact that three arbitrary matrices , and