Spin Diffusion Equations for Magnetized or Orbital Polarized Systems

TL;DR

This paper improves the derivation of spin diffusion equations in spintronics by using linear response theory and self-consistent approximations to ensure charge conservation in systems with spin-orbit interaction and polarization.

Contribution

It introduces a method combining linear response theory and self-consistent Born approximation to derive charge-conserving spin diffusion equations.

Findings

Standard coarse graining can violate charge conservation.

The proposed method satisfies the Ward-Takahashi identity.

Improved diffusion equations accurately model spin and charge transport.

Abstract

Charge and spin transport in spintronics devices can be described by a spin diffusion equation suitable for modelling scales much larger than the scattering and atomic scales. This work concerns the coarse graining procedure used to compute the coefficients of the diffusion equation, which are sensitive the details of individual atoms and impurities. We show with two simple examples that in spintronics devices which have both a spin-orbit interaction and magnetization, standard coarse graining can easily obtain diffusion equations which fail to conserve electronic charge. The same failure can occur in systems with both a spin-orbit interaction and orbital polarization. We show that linear response theory, coupled with the self-consistent Born approximation and ladder diagrams, offers an improved way of calculating diffusion equations. We show that the resulting equations satisfy a Ward-Takahashi identity that guarantees charge conservation.

Equations

$$\tau \frac{\partial }{\partial t}\vec{\rho }=D\vec{\rho }$$

$${\Sigma }_{\alpha \beta }^R(E)$$

$$\Sigma =-ı\hslash /2\tau$$

$$\hslash /2\pi \tau$$

$$I_{\gamma \zeta \lambda \delta }(\mathbf{q},E,\omega )$$

$${\overline{\Phi }}_{\alpha \nu \eta \beta }(\mathbf{q},\omega ,E)$$

$${\Phi }_{\alpha \nu \eta \beta }(\mathbf{k},\mathbf{q},\omega ,E)$$

$$\tau {\partial }_t\vec{\rho }$$

$$\vec{\rho }(t)$$

$$D_{\alpha \gamma \delta \beta }(\mathbf{q},\omega ,E)$$

$$((\hat{\chi }{)}^{-1}{)}_{\alpha \gamma \delta \beta }{\rho }_{\gamma \delta }$$

$${\hat{\chi }}_{\alpha \gamma \delta \beta }(\mathbf{q},\omega ,E)$$

$${\Sigma }_{\nu \eta }(E+\omega )-{\Sigma }_{\nu \eta }(E)$$

$${\hat{\Gamma }}_{\nu \zeta \lambda \eta }(\mathbf{q},\omega ,E)$$

$$H$$

$$n$$

$$\Gamma$$

$$\Sigma$$

$$G^R(\vec{k},s)$$

$$(D^{-1}(\vec{q}=0,\omega =0){)}_{ij}$$

$$(D^{-1}{)}_{00}$$

$${\Sigma }^A$$

$${\Sigma }^R$$

$$H^{A,R}$$

$$=\frac{k^2}{2m}+(\vec{b}+\vec{\Sigma })\cdot \vec{\sigma }+{\Delta }_{SO}k\left[\begin{array}{cc}0& e^{-ı{\varphi }_k}\\ e^{ı{\varphi }_k}& 0\end{array}\right]$$

$$=\frac{k^2}{2m}+(\vec{b}+\vec{\Sigma })\cdot \vec{\sigma }+{\Delta }_{SO}\left[\begin{array}{cc}0& (k_x-ık_y)\\ (k_x+ık_y)& 0\end{array}\right]$$

$$E^{A,R}(s,\vec{k})$$

$$E^{A,R,z}(s,\vec{k})$$

$$G^{A,R}(\vec{k},E)$$

$${\vec{G}}^{A,R}(\vec{k},E)$$

$${\vec{G}}^{A,R,z}(\vec{k},E)$$

$$X_1^{A,R}$$

$$X_2^{A,R}$$

$$X_3^{A,R}$$

$$\cos {\theta }_B^{A,R}$$

$$(E-(k^2/2m+{\Delta }_{SO}(k_x{\sigma }_x+k_y{\sigma }_y)+B_z{\sigma }_z){)}^{-1}=\frac{1}{2}\sum _s\frac{1+s\vec{X}(\vec{k})\cdot \vec{\sigma }}{E-E(s,\vec{k})}$$

Full text

I Introduction

Note: This paper is in a draft status, does not mention all relevant citations, and may still contain errors. We expect that any errors are relatively minor and do not affect the main message or results, and that all the major citations are already included.

I.1 What we want to model; the reason for this work

The immediate stimulus for this work was the recent fabrication of magnetic memory elements in which the magnetization switching is caused by an electric current, not a magnetic field. (See Miron et al, Nature 476, 189-193, and Buhrman et al, Science 336, 555-558.) Miron et al. (2011); Liu et al. (2012) These devices are bilayers, with a magnetic layer on top of metallic layer with a strong spin-orbit coupling. For instance, Miron’s device was a Cobalt-Platinum bilayer. When a electric current flows through these devices, the spin-orbit coupling causes spin currents and accumulations, and switches in the magnetization are observed. These devices offer possibilities both for new memory technologies and for integrating memory elements into CPUs. We emphasize that, unlike some other earlier spintronics devices, spin-orbit coupling and its effect on transport are key to the new devices.

Modeling these memory devices requires a microscopic description of the atoms in each of the layers, the spin-orbit coupling in the metallic layer, the interface between the metallic and ferromagnetic layers, and the magnetization in the magnetic layer. Impurity atoms, scattering, and mismatches between the layers are likely to be important. Along with this atomic scale physics, the modeling must be able to capture the much larger scale of the entire memory device - and the first devices were micron sized.

Similar modeling capabilities are required by many other candidates for new spintronics technologies: antiferromagnetic memories (Science v. 351 p. 587, 2016 Wadley et al. (2016)), graphene with adsorbed impurities producing a very large spin Hall effect (Nature Communications 5, 4748 Balakrishnan et al. (2014)), magnetic domain wall motion in racetrack memories, and devices involving the movement of magnetic skyrmions.

The results presented here concern the specific technical challenge of modeling devices which include both spin-orbit couplings and one of the following two properties: magnetization and/or orbital polarization.

I.2 Coarse Graining with Spin Diffusion Equations

Several modeling techniques have been developed to meet these modeling requirements. One approach is to build a model of the device which explicitly includes all the atoms in the device, and then to calculate its behavior using the Schodinger equation. This approach always involves a limit on the size of the device that can be simulated, but device sizes in excess of a billion atoms have been achieved using tight binding Hamiltonian representations of each atom - see for instance Ferreira’s KITE software (https://quantum-kite.com/resources/). Important work in the same direction has been done using ab initio (DFT) software, although the computational expense significantly reduces the system sizes which can be achieved. Some of the leaders in this direction, i.e. some of the people who have often focused specifically on ab initio calculation of spintronics devices and spin transport and dynamics, have been Paul Kelly, Ingrid Mertig, Stefan Blugel, the Prague research group of Kudrnovsky, Drchal, and Sob, and others.

Our work presented here concerns a different strategy for modeling spintronics devices, which is often called spin diffusion equations. Broadly speaking, this is a coarse-graining strategy. Starting with the microscopic description of the atoms within a sample, we determine a diffusion equation which controls charge and spin transport at distances much longer than the scattering length. The diffusion equation does describe the influence of atomic-scale processes on long-distance physics, but it does not include a description of transport at the atomic scale. In a model where each atom has only a spin up and a spin down state, and no additional states associated with different atomic orbitals, diffusion equations describe the evolution of four coarse grained variables: the charge density $\rho$ and the spin densities $S_{x},\,S_{y},\,S_{z}$. These are often written as a charge-spin density four-vector $\vec{\rho}=\left[\rho,\,S_{x},\,S_{y},\,S_{z}\right]^{T}$. In models with two or more atoms, or with several orbitals per atom, the charge-spin density vector has $4N^{2}$ components, where $N$ is the number of atomic orbitals per unit cell. The diffusion equation is then written as

[TABLE]

The diffusion operator $D$ is a $4N^{2}\times 4N^{2}$ operator, and is derived from the microscopic description of the model’s atoms and the disorder. It includes the physics of magnetic fields, spin precession, spin-charge conversion, and spin relaxation in both the Elliot-Yafet and the Dyakonov-Perel regimes. $\tau$ is the momentum relaxation time caused by scattering. Since this is a first-order partial differential equation, the initial value of $\vec{\rho}$ must be determined, and also geometry and boundary conditions of the sample, and then the spin diffusion equation can be solved straightforwardly.

Spin diffusion equations have four principal advantages: (1) For simple models the spin diffusion equation can be derived analytically, (2) the coefficients of the spin diffusion equation invite theoretical interpretation and analysis, (3) solution of the coarse-grained spin diffusion equation requires much fewer computational resources and therefore allows calculation of samples with realistic sizes, and (4) the coarse-grained charge and spin flow determined by a spin diffusion equation is easier to examine and comprehend than results describing spin and charge at each atomic site. There are also some disadvantages: (1) the spin diffusion equation approach involves approximations which typically include omitting certain quantum interference processes and also many-body electronic correlations, and (2) the process of starting from a microscopic description of a sample’s atoms and then deriving the correct spin diffusion equation requires the use of perturbation theory and Feynman diagrams, and therefore requires considerable care and expertise. The correct derivation of spin diffusion equations is the subject of ongoing research and investigation. In particular the new results presented here are about how to derive correct spin diffusion equations in materials which include both spin-orbit coupling and one of two following effects: magnetization and/or orbital polarization.

I.3 The Standard Techniques for Calculating the Spin Diffusion Equation

There are two well-known procedures for calculating the diffusion operator $D$:

• •

The Burkov - Nunez - Macdonald approach (PRB 70 155308 Burkov et al. (2004)), which is motivated by linear response theory.

• •

The Mishchenko - Shytov - Halperin approach (PRL 93 226602 Mishchenko et al. (2004)), which is noteworthy for starting within a Keldysh formalism.

While these two techniques have some significant differences, they have many points in common in their philosophy, development, approximations, and actual computational steps. Both approaches obtain identical results for some simple 2-D models with rotational in-plane symmetry, as was demonstrated already in the original publications.

The work presented in the present article was stimulated by the fact that both approaches (as originally presented) treat the momentum relaxation time $\tau$ as a single parameter. This is incorrect in systems with magnetization or orbital polarization, where different orbital or spin states can display different momentum relaxation times. Moreover, often $\tau$ is treated as an input parameter, without calculation from the microscopic quantities that are responsible for scattering. We will show that in materials with spin-orbit coupling and with also either magnetization or orbital polarization, this simplified treatment of $\tau$ is incorrect and produces incorrect results; $\tau$ must be calculated from the microscopic data instead of simply being chosen as an input paremeter. In our view neither the Burkov - Nunez - Macdonald approach nor the Mishchenko - Shytov - Halperin approach correctly handle the particular case of systems with both spin-orbit coupling and with also either magnetization or orbital polarization. We are not aware of work generalizing either of the two methods to handle this case. (In particular Wang and Manchon’s PRL 108 117201 on CoPt bilayers Wang and Manchon (2012) retained the assumption that $\Sigma$ is proportional to the identity.)

In the remainder of this article we will focus on the Burkov - Nunez - Macdonald approach (PRB 70 155308 Burkov et al. (2004)), which we will call ”the BNM approach”, although admittedly the literature is fairly evenly split between this approach and the Mishchenko - Shytov - Halperin approach Mishchenko et al. (2004). Since both approaches share so much in common, and since both display the same methodological failure when handling magnetized and/or orbital-polarized systems, our critique of and improvement on the BNM approach reflects also on the Mishchenko - Shytov - Halperin approach. In the linear response regime the Keldysh formalism, while more powerful for treating interacting systems out of equilibrium, cannot fail to give the same results as linear response theory, if both are calculated exactly. Therefore the only possible source of differences between the two approaches is from the approximations which each employ.

Here we briefly review the equations of the Burkov - Nunez - Macdonald approach , which determine both the self-energy $\Sigma$ and the diffusion operator $D$. We conflate the spin and orbital degrees of freedom into a single index which we denote by greek letters. We use the convention that the Green’s function is $G^{A,R}_{\alpha\beta}=(E-H\mp\Sigma^{A,R})^{-1}$, where $\Sigma^{A,R}$ is the self-energy. The $A,R$ superscript notation selects the advanced or retarded Green’s function $G$ and corresponding self-energy $\Sigma$. $\Sigma$ is subtracted for the advanced Green’s function and added for the retarded Green’s function. The weak Gaussian disorder in the potential $V$ is described by its second moment $\Gamma_{\alpha\gamma\delta\beta}\propto\langle V^{2}\rangle$, which has four indices because it is the second moment of $V$, which like $H$ has two indices. The self-consistent Born approximation (SCBA) which we use to approximate the self-energy is:

[TABLE]

Often this integral is approximated by keeping only the imaginary part of the Green’s function and treating it as a delta function confining the integration ${d\boldsymbol{k}}$ to the Fermi surface.

Both Burkov Nunez and Macdonald Burkov et al. (2004) and Mishchenko Shytov and Halperin Mishchenko et al. (2004) assume that $\Sigma_{\alpha\beta}$ is imaginary and proportional to the identity $1$. Since the imaginary part of $\pi^{-1}\int{d\boldsymbol{k}}G_{\gamma\delta}(\boldsymbol{k},E)$ is the single particle density of states, if the disorder has trivial spin structure i.e. $\Gamma\propto 1$, then the assumption that $\Sigma_{\alpha\beta}\propto 1$ is equivalent to assuming that the system has no spin or orbital polarization. Within this assumption, the self-energy is commonly written as

[TABLE]

where $\tau$ is the scattering time.

The SCBA equation is a self-consistent equation, since the Green’s function $G$ depends on the self-energy $\Sigma$. Solving this nonlinear equation requires iterative recalculation of $\Sigma$ from $G$, and $G$ from $\Sigma$, and as a rule must be done numerically rather than analytically. Therefore theoretical papers often avoid the numerical work and instead use

[TABLE]

It can be tempting to further simplify things by using the bare Green’s function $G(\Sigma=0)$ to evaluate $\rho$, thus avoiding the iterative self-consistency prescribed by the SCBA. Alternatively, while in principle remembering that $\rho$, $\Sigma$, and $1/\tau$ include self-consistency effects and all powers of $\Gamma$, some authors do not go through the steps of requiring numerical self-consistency, and leave these parameters as basically external parameters. In either case, the effect is to cut the connection between (A) the scattering processes that determine $\Sigma$ and the single-particle Green’s function and (B) the scattering processes that determine the diffusion operator $D$. When this connection is cut the resulting diffusion operator will fail to conserve charge if $\Sigma$ is not proportional to the identity, as we will show in the next sections. The only solution is to insist both analytically and numerically that $\Sigma$ be the real numerical solution of the SCBA for the specific values of $\Gamma$ and $G_{0}$ realized in the system one is modelling.

Burkov Nunez and Macdonald calculated a collision integral $I$ which represents one scattering event where both an electron and a hole scatter off the same impurity:

[TABLE]

Similarly to the integral determining the self-energy, the collision integral is often approximated by keeping only the imaginary part of one of the Green’s functions and treating it as a delta function confining the integration ${d\boldsymbol{k}}$ to the Fermi surface.

Burkov Nunez and Macdonald determine the diffusion operator $D$ directly from the collision integral:

[TABLE]

Actually one should expand $D$ to second order in momentum $\boldsymbol{q}$ and zeroth order in energy $\omega$ in order to obtain the diffusion operator. If one does not do so then it has a non-trivial time dependence and it contains all powers of the spatial derivatives. There isn’t much point in going to higher orders because the main benefit of going to higher order would be to investigate small length scales, but information about small length scales is removed and corrupted by BNM’s ladder approximation. Since this information is absent, attempts to recover it with higher derivatives produce results without any physical meaning.

The evolution operator $\exp(Dt/\tau)$ contains all powers of the collision integral $I$, representing repeated scatterings off of impurities. In perturbation theory the collision integral $I$ looks like one rung of a ladder, so powers of $I$ are called ladder diagrams, and the Burkov Nunez and Macdonald formula is a ladder approximation.

Some recent papers by Ado Dmietriev Ostrovsky and Titov Ado et al. (2015) and by Milletari and Ferreira Milletarì and Ferreira (2016) have shown that even though the ladder approximation calculates the longitudinal conductivity correctly, it is not sufficient for getting the transverse conductivity and spin Hall physics. The ladder approximation will give the correct order of magnitude for the Hall coefficient but can get the sign wrong, as well as the numerical coefficient. The solution to this problem is to renormalize the scattering tensor $\Gamma$ with a one loop integral.

I.4 Failure in Systems which are Magnetized or have Orbital Polarization

If the self-energy is a proportional to the identity, i.e. $\Sigma=\imath\hbar/2\tau$, then the BNM approach produces a diffusion operator $D$ which conserves the electronic charge component of $\vec{\rho}$. We will prove this statement in detail later in this article. The proof relies on $\Sigma$ commuting with the bare Hamiltonian $H_{0}$ inside the Green’s function, allowing $\tau$ it to be pulled outside of the integral over ${d\boldsymbol{k}}$. A similar trick can be played inside the momentum integral determining $\Sigma$. Using these tricks, one can show that the collision integral $I$ is independent of $\tau$ and therefore is able to cancel exactly against the $1$ in $D=1-I$ independently of the value of $\tau$. The cancellation provides an exact zero eigenvalue in $D$ which means that the inverse of the charge lifetime is zero; i.e. charge is conserved. Thus the commutation between $\Sigma$ and $H_{0}$ provides an elementary mathematical mechanism for guaranteeing charge conservation.

On the other hand, when $\Sigma$ is not proportional to identity, it in general does not commute with $H_{0}$, and its role in the collision integral $I$ and in $\Sigma$ becomes non-trivial. In this case the mechanism discussed above for guaranteeing charge conservation does not apply, and deeper physics/math must be used to obtain charge conservation. If a calculation of the diffusion operator $D$ omits or truncates the necessary physics, then the resulting $D$ will not conserve charge. We will show two very simple magnetized models which, when one uses the BNM approach and retains the fact that $\Sigma$ is not proportional to the identity, produce $D$’s which exhibit charge non-conservation. When the spin splittings in $H_{0}$ and in $\Sigma$ are comparable, charge is lost at a rate which is similar to the scattering rate $\tau^{-1}$. In this case the diffusion operator completely loses its validity.

I.5 The Correct Formalism

The solution to this difficulty is to ground the BNM approach more firmly in linear response theory. Rather than attempting to provide a general evolution equation describing the behavior of an arbitrary charge and spin distribution, we restrict our focus to a system’s response to small external perturbations $V$. Linear response theory allows one to compute $\rho$, the change in charge and spin density induced by $V$:

[TABLE]

Here we are using the Einstein convention where repeated indices have implicit sums. The key machinery in this equation is the susceptibility $\hat{\chi}$. The inverse susceptibility $(\hat{\chi})^{-1}$, if calculated exactly, would contain a complete description of the linear response at both short and long scales, including ballistic physics, diffusive motion, and Anderson localization. One can obtain from $(\hat{\chi})^{-1}$ the diffusion operator $D$ which governs long time and distance scales by expanding to second order in momentum $\boldsymbol{q}$ and first order in energy $\omega$.

In the ladder approximation, $(\hat{\chi})^{-1}$ has a simple relation to the diffusion operator $D=1-I$ calculated in the BNM approach:

[TABLE]

All of the details of calculating the collision integral $I$, $\overline{\Phi}$, and $\Sigma$ are the same as in the BNM approach. Charge conservation will not be recovered, however, unless the self-energy $\Sigma$ is the self-consistent solution of the SCBA equation.

The advantage of following linear response theory and calculating the susceptibility $\hat{\chi}_{0}$ is that a powerful Ward-Takahashi identity guarantees charge conservation. (todo: advanced-retarded, and check the explanation of $\Lambda$) The Ward-Takahashi identity reads:

[TABLE]

Here we are using Matsubara frequencies $\omega$ which are imaginary. Upon analytic continuation to the real axis the left hand side of the identity becomes $\Sigma^{A}_{\nu\eta}(E+\omega)-\Sigma^{R}_{\nu\eta}(E)$, which is in general non-zero when $\omega\rightarrow 0$. $\Lambda$ is the interacting vertex function. It can be obtained from susceptibility $\chi$ by pruning two Green’s functions and then substracting $1$, i.e. roughly $\chi=G^{A}G^{R}+G^{A}G^{R}\Lambda$. The $\omega$ multiplying the right-hand side of the Ward-Takahashi identity implies that at small $\omega$ (i.e. long times) $\Lambda$ and $\hat{\chi}$ both diverge like $\omega^{-1}$. This in turn implies that $(\hat{\chi})^{-1}$ has an eigenvalue which is zero when $\omega=0$. It is this eigenvalue which ensures that both $(\hat{\chi})^{-1}$ and the diffusion operator $D$ which can be obtained from it by using a Taylor series expansion conserve charge.

The Ward-Takahashi identity was proved in Ramazashvili’s PRB 66 220503(R), following similar earlier proofs Ramazashvili (2002). In this present work we generalize that proof to systems with magnetization or orbital polarization by keeping track of the spin/orbital indices at every step in the proof. The proof is fully rigorous and non-perturbative, requiring only that global charge be conserved. The price of this rigor is that one must restrict one’s scope to linear response theory, i.e. charge and spin response to small perturbing potentials.

Although the Ward-Takahashi identity guarantees that a system’s exact inverse susceptibility $(\hat{\chi})^{-1}$ conserves charge, it does not make guarantees about approximations to the susceptibility, including the ladder approximation $\hat{\chi}_{\alpha\gamma\delta\beta}=\sum_{\nu\eta}\overline{\Phi}_{\alpha\nu\eta\beta}((1-I)^{-1})_{\nu\gamma\delta\eta}$ which we have just given. Therefore the present work presents several checks that the ladder approximation, in combination with a full self-consistent SCBA calculation of the self-energy, does obey the Ward-Takahashi identity. Within the ladder approximation to $\hat{\chi}$, the Ward-Takahashi identity is

[TABLE]

I.6 What can we gain?

• •

Most importantly, we gain a working coarse-grained theory of spin transport that can handle modern spintronics devices much larger than the atomic scale.

• •

Use of the SCBA to perform a correct calculation of the self-energy results in a renormalization of the exchange field $b_{z}$ and the total magnetization. Moreover the Fermi velocity, which figures in the diffusion constants, spin-spin couplings, and charge-spin couplings, is renormalized.

• •

Comparing to Burkov Nunez and Macdonald’s diffusion operator $D=1-\Gamma\overline{\Phi}$, the new diffusion operator is $(1-\Gamma\overline{\Phi})\overline{\Phi}^{-1}$. In the regime Dyakanov-Perel spin relaxation, where the spin relaxation rates $\propto(1-\Gamma\overline{\Phi})\tau^{-1}\ll\tau^{-1}$, the new matrix $\overline{\Phi}^{-1}$ can safely be truncated at zeroth order in $\omega$ and $\boldsymbol{q}$. Therefore its only effect is to rotate the Burkov Nunez Macdonald kernel $(1-\Gamma\overline{\Phi})$, causing the spin-charge density $\vec{\rho}$ to be rotated. In the Elliot-Yafet regime this approximation is not valid and the new matrix must retain its $\omega$ and $\boldsymbol{q}$ dependence, causing changes in the spin diffusion operator (compared to the BNM operator) that go well beyond rotations of $\vec{\rho}$.

• •

The new introduced matrix $\overline{\Phi}^{-1}$ may be nontrivial even when there is no magnetization. Where this is so, past published results on spin diffusion equations in the Elliot-Yafet regime will need revision.

• •

Where the BNM approach uses a simple time derivative $\tau\frac{\partial}{\partial t}$, the correct treatment requires expanding $(\hat{\chi})^{-1}$ to first order in $\omega$, and $\frac{d(\hat{\chi})^{-1}}{d\omega}$ multiplies the time derivative, replacing $\tau$. Therefore the diffusion equation’s time derivative acquires a matrix structure which can include physics well beyond that of $\tau$. This non-trivial matrix structure of the time derivative can occur in non-magnetized systems, and if so would force revision of previous literature.

I.7 Outline of this Paper

The remainder of this work is organized as follows:

• •

Section II shows that the BNM approach produces diffusion operators which do not conserve charge when applied to two simple models: (a) a minimal model with magnetization and spin-orbit disorder, and (b) a model of CoPt which can also be applied to magnetized TIs.

• •

Section III contains (a) a general proof of the Ward-Takahashi identity, (b) derivation of the transport equation written in terms of $\hat{\chi}$, both exactly and when using the ladder diagram approximaiton, and (c) notes on how to perform the analytic continuation from Matsubara frequencies to real frequencies and energies. Much of the analytical work in this paper, and in particular the proof of the Ward-Takahashi identity, is done using Matsubara frequencies.

• •

Section IV verifies that the Ward-Takahashi identity is obeyed when using the SCBA plus ladder diagrams. It does so (a) in any model to first and second order in $\Gamma$, (b) to all orders when $\Sigma$ is proportional to the identity, (c) to all orders in a minimal model with magnetization and spin-orbit disorder - the same one which we showed does not conserve charge in the BNM approach, and (d) to second order in $\Gamma$ a model of CoPt which can also be applied to magnetized TIs - the same one which we showed does not conserve charge in the BNM approach. The derivations to second order can be extended to all orders by following the logic of Vollhardt and Wolfle’s PRB 22 4666 Vollhardt and Wölfle (1980).

• •

Section V has notes toward calculating the coefficients of the diffusion operator for a simple spin $1/2$ model with both magnetization and a spin-orbit interaction. This model has been used to analyze CoPt bilayers, and is also applicable to magnetized TIs. The work reported in this section is analytical, and is incomplete because it remains to solve the SCBA self-consistent equation numerically and then evaluate all the spin diffusion coefficients numerically.

• •

Section VI has some of our notes on the Bethe-Salpeter equation, Quantum Kinetic Equation, and charge conservation in this context. We worked on this topic a lot because of its connection with Mischenko Shytov and Halperin’s Keldysh approach to spin diffusion.

II Two Models Which Exhibit Charge Non-Conservation

In this section we will present two simple spin 1/2 models, and show that when treated using the BNM approach they produce spin diffusion equations which do not conserve charge.

II.1 Minimal Model: Magnetization Plus Spin-Orbit Disorder

This minimal model has two orbitals (spin up and spin down), a $\hat{z}$ Zeeman splitting between its energy levels which causes in-plane precession of the spin, a minimal dependence on $\vec{k}$ i.e. no spin-orbit coupling, a density of states which is different for spin up than for spin down (the asymmetry is controlled by $n_{z}$, the spin-$z$ component of the density of states), and both a spin-conserving disorder $\Gamma_{0}$ and a spin-orbit disorder $\Gamma_{x}\sigma_{x}\otimes\sigma_{x}$. $\Gamma_{0}$ is not necessary for obtaining the charge non conservation. What is necessary to produce this pathology is the presence of two non-commuting operators, the spin-orbit disorder $\Gamma_{x}\sigma_{x}\otimes\sigma_{x}$ and the self-energy which has a component proportional to $\sigma_{z}$.

We present the Hamiltonian, its eigenvalues, and the density of states $n$ (without disorder). $H$ is the Hamiltonian, $E$ is its eigenvalues, $s=\pm 1$ chooses among two eigenvalues, $n$ is the density of states, $\Gamma$ is the scattering vertex, $\Sigma$ is the self-energy, and $G^{R}$ is the retarded Green’s function.

[TABLE]

Note that we have calculated the self-energy in a way that is often used in the literature, i.e. without making sure that it is the self-consistent solution of the SCBA equation. We now follow the BNM approach for calculating the diffusion operator $D^{-1}$ at zero momentum and zero frequency. This object controls decay and mixing of spin and charge.

[TABLE]

The sign $t=\pm 1$ is $+1$ if the $\Gamma_{x}$ scattering comes from spin-orbit impurities, and $-1$ if it comes from local magnetizations - see Hikami Larkin and Nagaoka’s original paper Hikami et al. (1980).

If we remove the spin-orbit disorder $\Gamma_{x}$, then we get

[TABLE]

This clearly conserves charge. There is no coupling between charge and spin, i.e. the upper row and left column are filled with zeros. Moreover charge does not decay - the upper left corner of the matrix is zero.

On the other hand if we remove the ordinary disorder $\Gamma_{0}$ and retain only the spin-orbit term $\Gamma_{x}$, then we get two terms:

[TABLE]

The second term clearly does not conserve charge: the number in the upper left corner indicates that charge decays, and the numbers in the upper right and lower left corners indicate that charge mixes with the $z$ component of spin $S_{z}$. The charge non-conservation effect is biggest, of order $O(1)$, when $\Sigma_{1}=\pi n_{1}\Gamma_{z}$ is comparable to $E_{z}$, i.e. when the self-energy caused by the spin-orbit term is comparable to the level splitting caused by the magnetization. The $O(1)$ magnitude obtained in this case means that charge decays at the same time scale as the momentum relaxation time $\tau$.

Mathematica allows us to do anything in between limits of $\Gamma_{0}=0$ and $\Gamma_{x}=0$. At lowest order in $\Gamma_{x}/\Gamma_{0}$, we find

[TABLE]

This result shows that the BNM approach produces charge non-conservation even without a polarized density of states $n_{z}$, a polarized lifetime $\Sigma_{z}$, or spin-conserving disorder $\Gamma_{0}$. It is enough to have (a) a non-zero Zeeman splitting $E_{z}$ and (b) spin-non-conserving disorder $\Gamma_{x}$.

II.2 CoPt Model: Magnetization with Spin-Orbit Coupling

In the minimal model we used spin-orbit disorder but no spin-orbit interaction. Here we do the reverse, i.e. no spin-orbit disorder and instead a spin-orbit interaction. The spin non-conservation is caused in this case by the non-commutation of the magnetization with the spin-orbit term. This model is directly applicable to Cobalt-Platinum bilayers, and has been studied by Pesin and Macdonald’s PRB 86 014416 Pesin and MacDonald (2012) and by Wang and Manchon’s PRL 108 117201 Wang and Manchon (2012). However previous works did not notice or address the fact that the self-energy is spin dependent, and that correct treatment of the self-energy results is charge non-conservation.

This model can also be used to study topological insulators with magnetization, by setting the $k^{2}/2m$ term in the Hamiltonian and the eigenvalues to zero.

In the following we often write the magnetic field as $\vec{b}$ thus allowing it to point in any direction, but we also at times assume that it is a perpendicular magnetic field $\vec{b}=b_{z}\hat{z}$ because this is simpler. We include a spin-orbit interaction $\alpha(k_{y}\sigma_{x}-k_{x}\sigma_{y})$. Because of the magnetic field, the density of states $\hat{\rho}$ and the self-energy $\Sigma$ are not proportional to the identity. Significant spin-polarization occurs only near the bottom of the dispersion, where only one of the two spin-split states coincides with the Fermi energy.

We will show that that the BNM approach results in charge non-conservation and a coupling between charge and spin.

The self-energy is:

[TABLE]

We are violating convention here since we have not written the $\imath$ in $\Sigma$ explicitly. In general $\Sigma$ is expected to have an important imaginary part, and the real part is usually neglected. Also note that if $b_{x}=b_{y}=0$, then also $\Sigma_{x}=\Sigma_{y}=0$.

The Hamiltonian plus self-energy is:

[TABLE]

Since $\vec{b}$ always occurs in combination with $\vec{\Sigma}^{A,R}$, we write the combination as $\vec{B}^{A,R}=\vec{b}+\vec{\Sigma}^{A,R}$.

The energies are as follows. Note that if there is a self energy and it has an imaginary part, then these energies are, in principal, complex.

[TABLE]

The Green’s functions are:

[TABLE]

The above formula for the Green’s functions is based on the following inversion formula:

[TABLE]

we prove the inversion formula as follows:

[TABLE]

Now we introduce the scattering tensor, which describes nonmagnetic white noise disorder:

[TABLE]

We write the SCBA equation for calculating the self-energy:

[TABLE]

If $\vec{b}=b_{z}\hat{z}$ then it is trivial to prove that $\Sigma_{x}=\Sigma_{y}=0$, and the remaining integrals have rotational invariance resulting in one dimensional integrals:

[TABLE]

This SCBA equation should be solved self-consistently. However, many authors do not obtain self-consistency and instead do the following zeroth order approximation, using the bare Green’s function. We follow this practice to illustrate the consequences:

[TABLE]

We begin calculating the collision integral $I$:

[TABLE]

Since we want to show that charge is not conserved in the BNM approach, we choose to focus on $I_{ij}(E,\omega=0,q=0)$. $I_{00}$ determines the charge lifetime, and $I_{0j}$ determines the coupling between charge and spin. We also insert the Green’s functions:

[TABLE]

Following standard practice, we now approximate the denominator, taking the imaginary part, and treating it as a Dirac delta function.

[TABLE]

Now we require that the magnetic field be perpendicular to the plane, i.e. $\vec{b}=b_{z}\hat{z}$, and take advantage of the resulting rotational symmetry:

[TABLE]

Now we specialize to calculating the charge lifetime, i.e. $I_{00}$. If $I_{00}$ is not equal to $1$, then charge has a lifetime.

[TABLE]

Now we analyze the deviation of $\vec{X_{A}}\cdot\vec{X_{R}}$ from $1$, which is its value in a system without spin polarization. We also break the sum $\sum_{s,\acute{s}}$ into contributions with $s=\acute{s}$ and contributions with $s=-\acute{s}$. (Here we are assuming $\rho_{A}=\rho_{R}$ to simplify the algebra, and later also $\Sigma^{R}=-\Sigma^{A}$, but the same argument goes through without these assumptions.)

[TABLE]

$\delta XX$ is independent of $\phi$ and the angular integral is trivial. The first of the four terms in $I_{00}$ would be exactly $1$ in systems without spin polarization. Here it deviates from $1$ meaning that charge is not conserved, because

[TABLE]

Next we evaluate $\delta XX$ and show that it is not zero, implying that charge is not conserved.

[TABLE]

This expression’s deviation from zero controls charge non-conservation. Its numerator (inside the square root) is proportional to $\Sigma^{A}_{z}\Delta_{SO}k_{F}$. Overall the deviation from zero is of order $1$ when $\Sigma_{z}$ and $\Delta_{SO}|k(s,\phi)|$ are of the same magnitude, i.e. when the spin-orbit splitting $\Delta_{SO}$ is of the same magnitude as the magnetization energy $\Sigma_{z}$. In this case we find that the decay constant is of order $O(1)$, meaning that the charge decay time is of the same order as the momentum relaxation time.

In summary, charge is not conserved, i.e. $I_{00}$ is not $1$, both because $\pi\sum_{s}\frac{\gamma\rho(s)}{2\,Im(E_{A}(s))}$ is not one, and because $\delta XX$ is not zero.

Lastly, we consider the charge-spin coupling, which should be zero.

[TABLE]

Now we decompose $X_{A},X_{R}$ into their mean value and difference. In a system without spin polarization, the difference $\delta X$ is zero. Because of the accompanying factor of $s+\acute{s}$, the mean value $X_{+}$ contributes only when $s=\acute{s}$. Similarly, the difference contributes only when $s=-\acute{s}$.

[TABLE]

In the first term, the angular integral evaluates to zero when $j=1,2$, but not when $j=3$; the coupling between charge and $S_{z}$ spin is non-zero. In the second and third terms $\delta X_{j}$ also integrates to zero when $j=1,2$, but not when $j=3$. In all terms $X_{3}$ is independent of $\phi$ so the angular integral is trivial.

In summary, the $j=3$ terms in $I_{0j}$ are non-zero, which means that charge mixes with $S_{z}$, which again means that charge is not conserved.

III Ward-Takahashi Identity and Diffusion Equations

In this section we derive a Ward-Takahashi identity which ensures charge conservation. Our work is inspired by two papers:

• •

Ward-Takahashi identities for disordered metals and superconductors, by Revaz Ramazashvili, PRB 66 220503(R), 2002 Ramazashvili (2002).

• •

Ward-Takahashi identity for nonequilibrium Fermi systems, by Velicky Kalvova and Spicka, PRB 77 041201(R), 2008 Velický et al. (2008).

We prove the Ward-Takahashi identity using very strong and general non-perturbative arguments which are articulated in these two papers. The proof applies to any interacting theory which commutes with the global charge operator. The main value-added of our presentation is that we are explicitly emphasizing spin (and orbital) indices, and their role in the Ward-Takahashi identity, where the other works did not. Our development was important for understanding how charge conservation can be guaranteed in spin-polarized or orbital-polarized systems. The Ward-Takahashi identity which we will prove is:

[TABLE]

$\Sigma$ is the self-energy. $\Lambda$ is a four-point function, and is dimensionless. This Ward-Takahashi identity applies to several variants of $\Lambda$: the advanced-retarded variant $\Lambda^{AR}$, the retarded-retarded variant $\Lambda^{RR}$, and the advanced-advanced variant $\Lambda^{AA}$.

Ramazashvilli’s PRB 66 220503(R) Ramazashvili (2002) showed that the Ward-Takahashi identity for $\Lambda^{AR}$ is equivalent to a statement that $\Lambda^{AR}_{\alpha\beta,\gamma\gamma}(E+\omega,E,p,p)$ has a pole at $\omega=0$, i.e.

[TABLE]

Moreover, if we keep the $k$ momentum in $\Lambda^{AR}$, we will get something like

[TABLE]

This is recognizable as a diffuson pole. The absence of a constant term in this expression, i.e. the absence of a lifetime, indicates that charge is conserved. From this form written as a pole one can derive a diffusion equation governing transport. The pole guarantees that the diffusion equation will satisfy charge conservation.

The Ward-Takahashi identity is very powerful for three reasons: (1) it concerns $\Lambda$ which is the unique operator that controls charge and spin transport (in the linear response approximation), (2) $\omega$’s explicit role in the Ward-Takahashi identity makes the existence of the pole $\Lambda$ unavoidable, and (3) the Ward-Takahashi identity is rigorous for every system which conserves global charge.

There is another Ward-Takahashi identity which has been proved by Vollhardt and Wolfle’s PRB 22 4666 Vollhardt and Wölfle (1980):

[TABLE]

We will call this Ward II, in contrast to the more powerful one which we favor in this text. Here $U$ is the irreducible four point vertex, and has units of energy squared. Vollhardt and Wolfle’s Ward II identity is proved perturbatively, and does not contain an explicit factor of $\omega$. We have not found it directly useful for proving charge conservation in the general case of spin-polarized or orbital-polarized systems. However, Vollhardt and Wolfle’s PRB 22 4666 Appendix A is very useful for understanding that in disordered systems the a ladder-diagram approximation for $\Lambda$, combined with the CPA/SCBA approximation for $\Sigma$, does satisfy charge conservation Vollhardt and Wölfle (1980). The perturbative expansion and manipulation used by Vollhardt and Wolfle’s PRB 22 4666 Appendix A to prove their Ward II identity transfers exactly to showing that CPA/SCBA + the ladder-diagram approximation satisfies charge conservation.

After proving the Ward-Takahashi identity, we derive transport equations. First we discuss the general case of a charge-conserving system, and then we present the equations for disordered systems with the CPA/SCBA + the ladder-diagram approximation. Lastly there is an unfinished section discussing how these results, derived using Matsubara frequencies, should be analytically continued to real frequencies.

III.1 The Ward-Takahashi identity

We will derive the Ward-Takahashi identity for the vertex correction under the imaginary-time (Matsubara) formalism, with $\hbar=1$. The action $S[\psi^{\dagger},\psi]$ is uniquely obtained from the Hamiltonian $H(\tau)[\psi^{\dagger},\psi]$ as

[TABLE]

with the fermionic field variables (Grassmann numbers) $\psi_{\alpha}^{({\dagger})}(\boldsymbol{r},\tau)$. The Hamiltonian $H[\psi^{\dagger},\psi]$ should be normal-ordered so that it can be appropriately converted to the path-integral formalism. The integral over the imaginary time $\tau$ is taken from [math] to $\beta=1/T$. With this action, the partition function is given by the path integral

[TABLE]

where we suppress the source terms for the field variables.

The quantum and thermal average of the observable $\mathcal{A}$ is given by the path integral

[TABLE]

The two-point Green’s function for the fermion is given by

[TABLE]

For the spinless non-interacting fermions, the action can be written in general as

[TABLE]

with the Fourier transformed field variables

[TABLE]

where $\omega_{n}=(2\pi/\beta)(n+1/2)$ is the fermionic Matsubara frequency. Thus the non-interacting partition function $Z$ can be straightforwardly obtained,

[TABLE]

where we have used the definition that Grassmann integral is equivalent to the derivative, and that the Grassmann derivatives also anticommute each other.

Similarly, the averaging over $\psi\psi^{\dagger}$ is given by

[TABLE]

As a result, the non-interacting Green’s function can be written in a simple form,

[TABLE]

There is only one momentum index and one frequency index because we are exploiting momentum conservation and energy conservation. In real space the Green’s function is

[TABLE]

We now generalize from scalar variables to a multi-component system, which is similar. The action, the partition function, and the Green’s function are given by

[TABLE]

respectively.

As long as the original Hamiltonian $H[\psi^{\dagger},\psi]$ conserves the total particle number, the action $S[\psi^{\dagger},\psi]$ is invariant under the time-independent global U(1) gauge transformation

[TABLE]

The conservation of total charge

[TABLE]

is equivalent to this global gauge symmetry.

We can also define the time-dependent U(1) gauge transformation

[TABLE]

with a time-dependent scalar function $\theta(\tau)$. Under this transformation, the variation of the action $S[\psi^{\dagger},\psi]$ comes only from a “Berry phase” term:

[TABLE]

We now derive the Ward-Takahashi identity related to the U(1) gauge invariance, namely total charge conservation. Since the path integral $\int\mathcal{D}\psi^{\dagger}\mathcal{D}\psi$ is invariant under the local (time-dependent) U(1) gauge transformation, the following equality is justified:

[TABLE]

Comparing the first and the fourth lines and taking $\theta(\tau)$ infinitesimal, we obtain

[TABLE]

Using the form of $\delta_{\theta}S$ in Eq. (129), we obtain the Ward-Takahashi identity relating the two-point and four-point functions,

[TABLE]

Here we have introduced the shorthand notation $\mathbf{x}=(\boldsymbol{r},\tau)$ and $\mathbf{K}\equiv(\boldsymbol{k},\omega_{n})$.

Next we Fourier transform the Ward-Takahashi identity from position-time coordinates $\mathbf{x}=(\boldsymbol{r},\tau)$ to momentum-frequency coordinates $\mathbf{K}\equiv(\boldsymbol{k},\omega_{n})$. Taking the specific form

[TABLE]

with $\bar{\omega}_{m}=(2\pi/\beta)m$ the bosonic Matsubara frequency, the left-hand and right-hand sides of Eq. (137) read

[TABLE]

We have defined the vertex function $\Lambda$, a four point function, which we will discuss in more detail in the next section.

As a result, the Fourier transformed Ward-Takahashi identity is

[TABLE]

Multiplying $[G^{-1}(\mathbf{K}+\mathbf{W}/2)]_{\alpha^{\prime}\alpha}$ from the left and $[G^{-1}(\mathbf{K}-\mathbf{W}/2)]_{\beta\beta^{\prime}}$ from the right, we obtain

[TABLE]

Using the specific form

[TABLE]

[TABLE]

we finally obtain

[TABLE]

which relates the self energy $\Sigma$ and the vertex correction $\Lambda$. This is the final form of the Ward-Takahashi identity. In the second line we have written the companion identity which applies to the $\Lambda$’s outer indices.

Note that the Ward-Takahashi identity proves that $\Lambda$ contains a pole at $\bar{\omega}_{m}=0$. This is true because the left hand side of the Ward-Takahashi identity does not go to zero when $\bar{\omega}_{m}\rightarrow 0$ - there is a branch cut in $\Sigma$ at $\bar{\omega}_{m}=0$. Therefore the right side of the Ward-Takahashi identity must also be non-zero. And therefore some component of $\Lambda$ must scale as $(\bar{\omega}_{m})^{-1}$ to cancel the explicit factor of $\bar{\omega}_{m}$.

If the disorder is static, independent of time, then the energy components of $\mathbf{K},\mathbf{K}^{\prime\prime}$ must be equal: $\omega=\omega^{\prime\prime}$. Therefore the sum $\omega$ or $\omega^{\prime\prime}$ in the above expressions is extraneous, and we need only a sum over $\mathbf{k}$ or $\mathbf{k}^{\prime\prime}$. The Ward-Takahashi identity is now

[TABLE]

These equations also prove that the susceptibility tensor $\chi$ contains a pole at $\bar{\omega}_{m}=0$, since $\chi$ is linearly related to $\Lambda$. Please refer to the next section which defines $\chi^{(0)}$ and $\chi^{(vc)}$. With those definitions in place, we obtain

[TABLE]

A similar identity shows that that there is also a pole in $\sum_{\mathbf{K}}\chi^{(vc)}_{\alpha\beta\gamma\alpha}(\mathbf{K};\mathbf{W};\mathbf{K}^{\prime\prime})$, i.e. the Ward-Takahashi identity applies both to $\chi$’s internal indices and to its external indices. Therefore we can conclude that the inverse of $\chi$ conserves charge, and is (up to an inconsequential multiplicative constant) the charge transport operator.

Both ${\chi}$ and ${\Lambda}$ are fully symmetric under interchange of inner and outer indices and retain the same pole.

This concludes our proof of the Ward-Takahashi identity.

III.2 Preparation for the Transport Equation

We define a local spin-charge density $\rho_{\alpha\beta}(\mathbf{x},\tau)$ and a perturbing local potential $V_{\alpha\beta}(\mathbf{x},\tau)$. The charge transport equation governing $\rho$ is

[TABLE]

The transport kernel is $\hat{\chi}$, which is related to the four-point function, namely the susceptibility tensor, by

[TABLE]

$\hat{\chi}^{-1}\rho$ has the same units as $V$, and $\chi$ has the meaning of $d\rho/dV$.

The diffusion operator, i.e. the equation which controls transport at scales larger than the scattering length/time, can be extracted by expanding $(\hat{\chi})^{-1}$ to linear order in $\omega^{\prime}$ and second order in $q$. After doing this expansion, if one Fourier transforms from $\omega^{\prime}$ to $t$ then one obtains a simple exponential evolution controlled by the kernel of $\hat{\chi}$.

We now discuss the Fourier transform of $\hat{\chi}$ in some detail, to get a good understanding of how the transport equations work. We require translational invariance and energy conservation, which makes the description in the reciprocal space efficient. The Fourier transform of the transport function $\hat{\chi}$, and of any other two point function such as the Green’s function $G_{\alpha\beta}$, is given by

[TABLE]

The Fourier transform of the susceptibility tensor, and of any other four point function, is

[TABLE]

We now separate the response function into two parts. The first part $\chi^{(0)}$ gives the complete response function of a non-interacting theory, and can be expressed as a product of two Green’s functions using Wick’s theorem. In interacting systems there is a vertex correction $\chi^{(vc)}$. We write $\chi^{(vc)}$ in terms of the vertex function $\Lambda$, a four point function. This is a definition of $\Lambda$:

[TABLE]

The Fourier transforms of these quantities are

[TABLE]

The Ward-Takahashi identity is

[TABLE]

where the self-energy is defined by

[TABLE]

We have neglected the diamagnetic contribution.

[TABLE]

III.3 The Diffusive Transport Equations - the Ladder diagram susceptibility and the SCBA

Now we introduce the ladder approximation for $\chi,\Lambda$. We also introduce the SCBA for $\Sigma$, which requires an assumption that disorder is weak and governed by Gaussian statistics. This is combination of approximations is exactly the same as the RPA with a local $\phi^{4}(x)$ interaction and an interaction strength set by $\Gamma$, with a single modification: the disorder is static, so most of the integrations over energy $\omega$ are removed. As long as the impurity density is small enough, we can easily treat strong impurities and non-Gaussian statistics by (a) moving from the SCBA to the CPA, and (b) adding T-matrices to the ladder approximation. See work by Aires Ferreira and collaborators for development of T-matrices Ferreira et al. (2014); Pachoud et al. (2014).

Within these approximations, transport equations are:

[TABLE]

$I$ is the collision integral. The diffusion operator can be extracted by expanding $(\hat{\chi})^{-1}$ to linear order in $\omega^{\prime}$ and second order in $q$. After doing this expansion, if one Fourier transforms from $\omega^{\prime}$ to $t$ then one obtains a simple exponential evolution controlled by the kernel of $\hat{\chi}$.

The Ward-Takahashi identity is

[TABLE]

Comparing to Burkov Nunez and Macdonald Burkov et al. (2004), their transport kernel $((\hat{\chi}^{transport})^{-1})$ is $1-I(\boldsymbol{q},E,\omega)$, where $I$ is a collision operator. In other words, they have, roughly, $\hat{\chi}^{transport}=(1-I(\boldsymbol{q},E,\omega))^{-1}$. The biggest differences between their result and ours are that (A) we have an additional factor of $\Gamma GG$ multiplying their kernel, and (B) we require that the equation for the self-energy $\Sigma$ be solved self-consistently.

III.4 Some notes on the Matsubara formalism

This section is not complete yet.

The goal of this section is to perform the transition from Matsubara frequencies to real frequencies. The discussion mostly follows the discussion in “Condensed Matter Field Theory” by Altland and Simons Altland and Simons (2010).

If we want to analytically continue to the real axis then the $\imath\omega$ prefactor turns into $\omega$, the quantities at $E$ become retarded quantities, and the quantities at $E+\omega$ become advanced quantities. One obtains:

[TABLE]

Suppose the Hamiltonian for the electrons is given by the matrix form,

[TABLE]

Then the imaginary-time action is given by

[TABLE]

where $\beta=(k_{B}T)^{-1}$ and $\omega_{n}=(2\pi/\beta)(n+1/2)$. We also add a source term coupled to generalized density,

[TABLE]

where $\bar{\omega}_{m}=(2\pi/\beta)m$ and

[TABLE]

The overall partition function is defined by the path integral,

[TABLE]

The expectation value of the charge is given by differentiating the partition function by $\phi$,

[TABLE]

Thus the response function $\chi$ is given by the density-density correlation,

[TABLE]

In the Fourier space, we have the familiar structure from the Wick’s theorem,

[TABLE]

The disorder-averaged Green’s function is given by

[TABLE]

As $h$ and $\gamma$ are Hermitian matrices, the Green’s function can be diagonalized by a certain unitary transformation $U$,

[TABLE]

The summation over Matsubara frequencies can be transformed into the loop integral in the complex plane using the Fermi-Dirac function $f(z)=(1+\exp(\beta z))^{-1}$,

[TABLE]

The contour must contain all the Matsubara frequencies $\imath\omega_{n}$ but it may not contain the branch cuts at $\mathrm{Im}z=0$ and $\mathrm{Im}(z+i\bar{\omega}_{m})=0$. Moreover parts of the contour at infinity are zero. The remaining integrals sum to

[TABLE]

By the analytic continuation $i\bar{\omega}_{m}\rightarrow\omega$, and then reorganizing, we obtain

[TABLE]

At this point one can take the imaginary part of the $\chi$ and after a few short steps obtain the Kubo formula for the conductivity, as shown in Bruus and Flensberg Bruus and Flensberg (2004).

IV Checks that our formalism combining the SCBA with a susceptibility built from ladder diagrams obeys the Ward-Takahashi identity

Although any field theory that conserves global charge must obey the Ward-Takahashi identity and charge conservation, approximations made to that theory may result in charge non-conservation if the approximation used to calculate the self energy $\Sigma$ does not match the approximation used to calculate the susceptibility $\chi$. We have proposed to use the SCBA for the self energy and a sum of ladder diagrams (very similar to RPA) for the susceptibility. In this section we perform many checks that the proposed approximations do satisfy the Ward-Takahashi identity.

Within the combination of SCBA and ladder diagrams, the Ward-Takahashi identity is:

[TABLE]

Sums are implied when we repeat indices.

In the next subsections discuss (A) checks of the Ward-Takahashi identity for all hamiltonians $H$ and disorder $\Gamma$, order by order in $\Gamma$. Next we show (B) that when the self-energy $\Sigma$ is proportional to the identity, the Ward-Takahashi identity is always satisfied. Next we show that (C) the minimal model (magnetization plus spin-orbit disorder) always satisfies the Ward-Takahashi identity. Lastly we develop (D) the CoPt model (magnetization with spin-orbit coupling) and show how to check the Ward-Takahashi identity order by order. In an earlier section we have showed that both the minimal model and the CoPt model result in equations which violate charge conservation when treated with the BNM approach, so the results shown here confirm that SCBA + ladder diagrams is superior.

IV.1 Check of the Ward-Takahashi identity at first and second order in the disorder $\Gamma$

At first order in $\Gamma$ the Ward-Takahashi identity is:

[TABLE]

So the Ward-Takahashi identity is satisfied at first order.

At second order in $\Gamma$ the Ward-Takahashi identity is:

[TABLE]

So the Ward-Takahashi identity is satisfied at second order.

@e believe that similar results can be proved at any order if one uses the strategy described by Vollhardt and Wolfle’s PRB 22 4666 Appendix A Vollhardt and Wölfle (1980). The examples we showed for first and second order already implement that strategy.

IV.2 Charge conservation when $\Sigma$ is proportional to the identity

Now we look at the simplifying case where $\Sigma$ is proportional to the identity, and prove that probability conservation is obtained trivially. In this case:

[TABLE]

The last line proves that the operator $1-I$, the diffusion kernel proposed by Burkov Nunez and Macdonald, conserves charge when $\Sigma$ is proportional to the identity. However it only proves charge conservation when the charge operator acts on one side of $I$, not when it acts on the other side of $I$. We can expect that if $\Gamma$ has spin structure then as a general rule $I$ will conserve charge on only one side.

In the last line we have assumed that $0\neq\lim_{\omega\rightarrow 0}-\Sigma(E)+\Sigma(E+\omega)$, which is reasonable because $\Sigma(E)$ is analytically continued to $\Sigma^{R}$ and $\Sigma(E+\omega)$ is analytically continued to $\Sigma^{A}$.

Continuing through with the full Ward-Takahashi identity in the case where $\Sigma$ is proportional to the identity:

[TABLE]

Q.E.D.

IV.3 Minimal Model: Magnetization Plus Spin-Orbit Disorder

We already discussed this model and showed that when treated with the BNM approach the resulting spin diffusion equations do not conserve charge. Here we show that using the SCBA plus ladder diagrams in the correct way one does obtain charge conservation.

In the minimal model there are two orbitals (spin up and spin down), $H$ is the Hamiltonian, $n$ is the density of states, $\Gamma$ is the scattering vertex, $\Sigma$ is the self-energy. This minimal model has a $\hat{z}$ Zeeman splitting between its energy levels which causes in-plane precession of the spin, has minimal dependence on $\vec{k}$, has a density of states which is larger for spin up than for spin down (the asymmetry is controlled by $n_{z}$), and has a spin-conserving disorder $\Gamma_{0}$ and a spin-orbit disorder $\Gamma_{x}\sigma_{x}\otimes\sigma_{x}$. $\Gamma_{0}$ is not necessary for obtaining the probability non conservation.

[TABLE]

We use the SCBA/CPA equation to calculate $\Sigma$:

[TABLE]

For ordinary disorder, the collision integral is:

[TABLE]

The $\sigma_{x}$ disorder just adds a $\sigma_{x}$ after (or before?) each of the Green’s functions, so for our model we use the equations below. At this point we introduce $t$, which is an extra sign found in the HLN paper. $t=+1$ if the $\Gamma_{x}$ scattering comes from spin-orbit impurities, and $t=-1$ if it comes from local magnetizations - see Hikami-Larkin-Nagaoka’s original paper Hikami et al. (1980).

[TABLE]

To work toward the Ward-Takahashi identity, we set $q=0$:

[TABLE]

Analyzing the matrix elements in the numerator gives:

[TABLE]

Converting to a different basis the identity and $\sigma_{z}$:

[TABLE]

There is also a $S_{x},S_{y}$ sector, but this sector is decoupled from the $1,S_{z}$ sector and doesn’t affect probability conservation.

In the $1,S_{z}$ sector we have $s=\acute{s}$, and $I$ simplifies to:

[TABLE]

In summary (and we can absorb the $(-2\pi m)$ factors into $\Gamma_{1},\Gamma_{z}$, giving $\Gamma$ units of energy):

[TABLE]

The susceptibilities are:

[TABLE]

Note that $1-I$ does not conserve charge, i.e. it is not zero when multiplied by either $\begin{bmatrix}1&1\end{bmatrix}$ (on the left) or $\begin{bmatrix}1\\ 1\end{bmatrix}$ (on the right). This shows again that the BNM approach for calculating spin diffusion equations does not satisfy charge conservation.

The Ward-Takahashi identity is:

[TABLE]

We check the Ward-Takahashi identity when $\Gamma_{z}=0$:

[TABLE]

So the Ward-Takahashi identity is satisfied when $\Gamma_{z}=0$.

We check the Ward-Takahashi identity when $\Gamma_{1}=0$:

[TABLE]

So the Ward-Takahashi identity is satisfied when $\Gamma_{1}=0$.

Moreover, using mathematica, we have checked that the Ward-Takahashi identity works for all values of $\Gamma_{1},\Gamma_{z}$. So the Ward-Takahashi identity is satisfied for all $\Gamma_{1},\Gamma_{z}$.

Investigating $1-I^{-1}$, we find that at $\omega=0$ one of its eigenvalues is zero. The matrix is not symmetric. The right eigenvector corresponding to the zero eigenvalue has a simple form.

We think that a spin diffusion equation does not have to be symmetric, but perhaps its $q=0,\omega=0$ component, i.e. its lifetime matrix, should be symmetric.

IV.4 CoPt Model: Magnetization with Spin-Orbit Coupling

This model can also be used to study topological insulators with magnetization, by setting the $k^{2}/2m$ term in the Hamiltonian and the eigenvalues to zero.

In the following we sometimes work with the magnetic field $\vec{b}$ pointing in any direction, but we also sometimes specialize to perpendicular magnetic field $\vec{b}=b_{z}\hat{z}$ because this is simpler. This model includes a spin-orbit interaction $\alpha(k_{y}\sigma_{x}-k_{x}\sigma_{y})$. The density of states $\hat{\rho}$ and the self-energy $\Sigma$ are not proportional to the identity. Significant spin-polarization occurs only near the bottom of the dispersion, where only one of the two spin-split states coincides with the Fermi energy.

The self-energy is:

[TABLE]

We are violating convention here since we have not written the $\imath$ in $\Sigma$ explicitly. In general $\Sigma$ is expected to have an important imaginary part, and the real part is usually neglected. Also note that if $b_{x}=b_{y}=0$, then also $\Sigma_{x}=\Sigma_{y}=0$.

The Hamiltonian plus self-energy is:

[TABLE]

Since $\vec{b}$ always occurs in combination with $\vec{\Sigma}^{A,R}$, we write the combination as $\vec{B}^{A,R}=\vec{b}+\vec{\Sigma}^{A,R}$.

The energies are as follows. Note that if there is a self energy and it has an imaginary part, then these energies are, in principal, complex.

[TABLE]

The Green’s functions are:

[TABLE]

The above formula for the Green’s functions is based on the following inversion formula:

[TABLE]

Now we introduce the scattering tensor, which describes nonmagnetic white noise disorder:

[TABLE]

We compute the self-energy with the SCBA/CPA:

[TABLE]

If $\vec{b}=b_{z}\hat{z}$ then it is trivial to prove that $\Sigma_{x}=\Sigma_{y}=0$, and the remaining integrals have rotational invariance resulting in one dimensional integrals:

[TABLE]

We begin calculating the collision integral $I$:

[TABLE]

We now resolve $I$ into its $1,S_{x},S_{y},S_{z}$ components by tracing over Pauli matrices. After doing the angular integral, the $S_{x},S_{y}$ components decouple from the $1,S_{z}$ components. We don’t evalulate the $S_{x},S_{y}$ components here because they do not affect the Ward-Takahashi identity.

Resolved into $1,S_{z}$ components, after performing the angular average, the numerator of $I$ is:

[TABLE]

When $\Sigma=0$, we obtain $\theta_{B}^{A}=\theta_{B}^{R}$, and the numerator simplifies enormously. When $s=\acute{s}$, the numerator is

[TABLE]

When $s=-\acute{s}$, the numerator is

[TABLE]

The fact that there is no coupling from charge to spin when $s=-\acute{s}$ is crucial for obtaining the Ward-Takahashi identity.

[TABLE]

The Ward-Takahashi Identity is:

[TABLE]

These equations are ready to evaluate numerically. The integrals do not have to be done right, as long as the integral in the self-energy matches the integral in $I$. More precisely, when $s=\acute{s}$ in $I$ and $\Sigma=0$ in $I$, then $-\omega I$ should be equal to the zeroth-order approximation to $\Sigma^{A}(E)-\Sigma^{R}(E)$. In particular, probably one could make a Fermi-surface approximation, though some care would have to be given to choosing $\tau$ for each value of $s$, and there would be some ambiguity about interpreting this choice. There is no numerical reason to make this approximation, it’s not any trouble to perform the integration, but in terms of physical interpretation it would allow some interpretation of the total result in terms of a density of states and a lifetime.

It is required that that the self-energy calculation be taken all the way to self-consistency, and that the calculation of $I$ use the correct self-energy. We also think that this means that it is not OK to approximate the difference in energies in the two denominators as $\omega+\imath/\tau$ - you have to use the correct self energies, and the correct $s,\acute{s}$.

At leading order in $\Gamma$ the Ward-Takahashi Identity is:

[TABLE]

At second order in $\gamma$ the Ward-Takahashi identity is again satisfied:

[TABLE]

We believe that similar results can be proved at any order if one uses the strategy described by Vollhardt and Wolfle’s PRB 22 4666 Appendix A Vollhardt and Wölfle (1980). The examples we showed for first and second order already implement that strategy.

V Spin Diffusion Equations for the CoPt Model: Magnetization with Spin-Orbit Coupling

This section has notes toward calculating the coefficients of the diffusion operator for a simple spin $1/2$ model with both magnetization and a spin-orbit interaction. This model has been used to analyze CoPt bilayers, and is also applicable to magnetized TIs by setting the $k^{2}/2m$ term in the Hamiltonian and the eigenvalues to zero. The work reported in this section is analytical, and is incomplete because it remains to solve the SCBA self-consistent equation numerically and then evalulate all the spin diffusion coefficients numerically.

In the following we work with $\vec{b}$ pointing in any direction, but we also keep the equations for perpendicular magnetic field $\vec{b}=b_{z}\hat{z}$ because this is simpler. We include a spin-orbit interaction $\alpha(k_{y}\sigma_{x}-k_{x}\sigma_{y})$. The density of states $\hat{\rho}$ and the self-energy $\Sigma$ are not proportional to the identity.

The self-energy is:

[TABLE]

We are violating convention here since we have not written the $\imath$ in $\Sigma$ explicitly. In general $\Sigma$ is expected to have an important imaginary part, and the real part is usually neglected. Also note that if $b_{x}=b_{y}=0$, then also $\Sigma_{x}=\Sigma_{y}=0$.

The Hamiltonian minus self-energy is:

[TABLE]

In the case of a TI we set $H_{0}(\boldsymbol{k}^{2})=0$. In the case of a 2-D electron gas we set $H_{0}(\boldsymbol{k}^{2})=h_{0}\boldsymbol{k}^{2}/2k_{0}^{2}$.

The eigenvalues are as follows. Note that if there is a self energy and it has an imaginary part, then these eigenvalues are, in principal, complex. The first equation gives the general case for any $\vec{b},\Sigma$, and the second equation gives the eigenvalues when $\vec{b},\Sigma$ are aligned with the $z$ axis.

[TABLE]

In the following expressions for the Green’s functions, we use the following inversion formula:

[TABLE]

The Green’s functions are:

[TABLE]

The approximation in the last line is justified because if we write $\boldsymbol{k}+\boldsymbol{q}/2$, or $E+\omega/2$, these have order of magnitude $q/k_{F}=l_{F}/L$, or $E/E_{F}=t_{F}/T$, which are small in the diffusive regime. It’s only when these are divided by $\Sigma/E_{F}=\hbar/(E_{F}\tau)=1/(lk_{F})=l_{F}/l\ll 1$, which are small if disorder is ”weak”, that we obtain quantities that are significant at diffusive length and time scales. The only way to get this extra multiplicative factor is from the denominator of the Green’s functions. Moreover only terms in the denominator that are linear in $\boldsymbol{q},\omega$ need to be considered, again because you can get only one compensating factor of $\Sigma/E_{F}$ from the denominator. In summary, the numerator’s dependence on $\boldsymbol{q},\omega$ should be omitted, and the denominator should be expanded only to linear order in these quantities.

We compute the self-energy with the SCBA:

[TABLE]

This equation must be solved iteratively until self consistency is obtained; otherwise the Ward-Takahashi identity will be broken.

If $\vec{b}=b_{z}\hat{z}$ then it is trivial to prove that $\Sigma_{x}=\Sigma_{y}=0$, and the remaining integrals have rotational invariance resulting in one dimensional integrals:

[TABLE]

We compute the main integral, assuming that the disorder has no spin structure:

[TABLE]

At this point we can do all the angular integrals. Writing $\overline{\Phi}$ now in terms of its spin and charge components, we obtain:

[TABLE]

In the next equations we redo the calculations with disorder that is purely magnetic, of the type where $M$ is proportional to $G^{A}\sigma_{x}\otimes G^{R}\sigma_{x}+G^{A}\sigma_{y}\otimes G^{R}\sigma_{y}$:

[TABLE]

It is possible to mix the disorder types, in which case $M$ will be a linear interpolation between the results with magnetic disorder and the results with spin-independent disorder.

There are some exotic terms in $\overline{\Phi}$ caused by $\Sigma^{A}_{z}\neq\Sigma^{R}_{z}$:

• •

The charge component $\overline{\Phi}_{00}$ and the charge-$S_{z}$ coupling $\overline{\Phi}_{03}=\overline{\Phi}_{30}$ acquire a part which is not proportional to $1+s\acute{s}$, allowing $s=-\acute{s}$ physics to contribute.

• •

The energy beating factor $\delta$, its derivative $W$, and the Fermi velocity $V$ acquire a part which is proportional to $s+\acute{s}$, allowing $s=\acute{s}$ physics to couple to the spin-orbit interaction.

• •

The charge-$S_{x},S_{y}$ couplings become nonsymmetric, i.e. $\overline{\Phi}_{01}-\overline{\Phi}_{10}$, $\overline{\Phi}_{02}-\overline{\Phi}_{20}$ become nonzero.

• •

The symmetric part of the charge-$S_{x},S_{y}$ couplings, i.e. $\overline{\Phi}_{01}+\overline{\Phi}_{10}$ and $\overline{\Phi}_{02}+\overline{\Phi}_{20}$, acquire a part which is not proportional to $s+\acute{s}$, allowing $s=-\acute{s}$ physics to contribute.

• •

The spin-mixing terms $\overline{\Phi}_{12}-\overline{\Phi}_{21}$, $\overline{\Phi}_{23}-\overline{\Phi}_{32}$, $\overline{\Phi}_{13}-\overline{\Phi}_{31}$ acquire a part which is not proportional to $s-\acute{s}$, allowing $s=\acute{s}$ physics to contribute.

VI Notes on the Bethe-Salpeter Equation and the Quantum Kinetic Equation

The Bethe-Salpeter equation is exact and offers a Ward-Takahashi identity, so we explored it extensively to see how it could help us with magnetized systems. One of the most attractive things about this equation is that it offers a nice base for developing and understanding the Keldysh based approach of Mischenko Shytov and Halperin, which is closely tied to a Quantum Kinetic Equation. Therefore after exploring the Bethe-Salpeter equation and its Ward-Takahashi identity we spent a lot of time trying to understand the Keldysh based approach and how it could be made fully general including for magnetized systems. Pletyukhov’s PRB 75 155335 was particularly helpful here Pletyukhov (2007). We were never able to bring this direction to a conclusion, and we still do not understand how to correctly generalize the Keldysh approach to magnetized systems.

Eventually we decided that the Ward-Takahashi Identity offered by the Bethe-Salpeter equation is not powerful enough, after we realized that the Ward-Takahashi identity explained by RamazashviliRamazashvili (2002) is more powerful. We did not figure out how to conserve probability in a magnetized system until we had abandoned the Bethe-Salpeter approach in favor of Ramazashvili’s approach.

Probably the reason why the Bethe-Salpeter equation approach did not give us what we need is that it naturally results in a description of how the density matrix evolves; a quantum kinetic equation whose inputs and outputs are both densities. This turns out to be a bit too ambitious. What we have succeeded in doing is producing a formalism which, given a perturbing potential as an input, produces a response density. We leave the problem of calculating the perturbing potential to another day. This approach allows us to build on the foundation of linear response theory, which is rock solid and gives us what we need to obtain spin conservation.

It is clear that there is a unique way within linear response theory (and still unique after application of the SCBA) to obtain charge conservation in magnetized systems. Therefore a correct generalization of the Keldysh approach must necessarily give identical results to the approach we present in this paper.

We retain these notes because the Bethe-Salpeter equation and the Keldysh approach retain great interest in their own right.

VI.1 Density matrix and probability conservation

Let us take an $N$-band system. With an arbitrary basis $|\mu\boldsymbol{k}\rangle\ (\mu=1,\ldots,N)$, the (unperturbed) Hamiltonian can be written in the $N\times N$ matrix form, $\hat{H}(\boldsymbol{k})=H_{\mu\nu}(\boldsymbol{k})$. Incorporating the disorder term $\hat{V}$, the time evolution of the density matrix $\hat{\rho}(t)$ is given by

[TABLE]

where $\hat{\rho}_{0}$ is the density matrix at the initial time $(t=0)$. The matrix element becomes

[TABLE]

Averaring over the disorder, we obtain

[TABLE]

We need to calculate the tensor $\check{\Phi}(\boldsymbol{k},\boldsymbol{k}^{\prime},\boldsymbol{q};E,\epsilon)$.

The total probability at time $t$ is given by

[TABLE]

Using Eq.(443), we obtain

[TABLE]

where

[TABLE]

Thus, the probability conservation law is expressed by

[TABLE]

We need to show the last equality, i.e. $\epsilon P(\epsilon)$ is independent of $\epsilon$.

VI.2 Bethe-Salpeter equation and kinetic equation

The disorder-averaged 4-point propagator $\Phi^{\mu\nu}_{\mu^{\prime}\nu^{\prime}}$ satisfies the Bethe-Salpeter equation,

[TABLE]

where $U^{\alpha\beta}_{\alpha^{\prime}\beta^{\prime}}(\boldsymbol{k},\boldsymbol{k}^{\prime\prime},E,\epsilon)$ denotes the “irreducible” vertex function, and the bare propagator $\Phi_{0}$ is defined by

[TABLE]

In the matrix form, this equation can be simplified as

[TABLE]

We should be careful that the advanced Green’s function is transposed.

Let us now fractionalize the bare propagator $\check{\Phi_{0}}$ as in the single-band case. In the matrix form, it becomes

[TABLE]

where

[TABLE]

Multiplying $\check{\mathcal{E}}(\boldsymbol{k},E,\epsilon)$ to Eq.(450) from the left, we obtain the kinetic equation

[TABLE]

leading to

[TABLE]

Here, the left hand side is the form appearing in Eq.(447).

VI.3 Ward-Takahashi identity and probability conservation

Let us now substitute the right hand side of Eq.(458) to Eq.(446), and manipulate their contributions to $\epsilon P(\epsilon)$ one by one.

(A) $\delta\check{H}\check{\Phi}$ term:

[TABLE]

(B) $\delta\check{\Sigma}\check{\Phi}$ term:

[TABLE]

In the last line we have used the Ward-Takahashi identity

[TABLE]

(C) $\check{\mathcal{G}}$ term:

[TABLE]

In the last line we have shifted the variable $E$ by $\pm\epsilon/2$ respectively. Using

[TABLE]

we obtain

[TABLE]

which is independent of $\epsilon$.

(D) $\check{\mathcal{G}}\check{U}\check{\Phi}$ term:

[TABLE]

where we have used the time-reversal symmetry of the vertex function,

[TABLE]

Now we can see that $\epsilon P^{\mathrm{(D)}}(\epsilon)$ cancels with $\epsilon P^{\mathrm{(B)}}(\epsilon)$.

Therefore, we obtain

[TABLE]

leading to

[TABLE]

Acknowledgements.

We are greatly indebted to Y. Araki and A. H. Macdonald for their support throughout this work. We gratefully acknowledge stimulating discussions with S. Kettemann, P. Wenk, X. Dai, B. A. Bernevig, A. Petrovic, A. Ferreira, H. W. Kim, A. Ho, M. Titov, P. Ostrovsky, J. Wunderlich, T. Ziman, M. Miron, M. Chshiev, B. Dieny, S. Parkin, J. Katine, P. Mavropoulous, K. Flensberg, S. Blugel, P. Kelly, M. Sob, and D. Edwards. We thank POSTECH, the University of York, the Institut Laue-Langevin, and the Forschungzentrum Juelich for hospitality. Conversations with C. Mueller and C. Miniatura were very influential early in the project. We acknowledge support from EPSRC grant EP/M011038/1.