Iterative path-integral summations for the tunneling magnetoresistance in interacting quantum-dot spin valves
S. Mundinar P. Stegmann J. König S. Weiss







Theoretische Physik, Universität Duisburg-Essen and CENIDE, 47048 Duisburg, Germany




We report on the importance of resonant-tunneling processes on quantum transport through interacting quantum-dot spin valves. 
To include Coulomb interaction in the calculation of the tunneling magnetoresistance (TMR), we reformulate and generalize the recently-developed, numerically-exact method of iterative summation of path integrals (ISPI) to account for spin-dependent tunneling. 
The ISPI scheme allows us to investigate weak to intermediate Coulomb interaction in a wide range of gate and bias voltage and down to temperatures at which a perturbative treatment of tunneling severely fails. 


Iterative path-integral summations for the tunneling magnetoresistance in interacting quantum-dot spin valves
    S. Weiss
    December 30, 2023
=============================================================================================================





§ INTRODUCTION

Spintronic devices such as spin valves rely on magnetoresistive effects, in which the charge current depends on the magnetic configuration of the involved magnetic components.
Examples are the giant- (GMR) and tunnel- (TMR) magnetoresistance effect.
The TMR effect, first  observed by Julliere<cit.> as a spin-dependent resistance (or conductance) in magnetic tunnel structures, is directly related to the degree of spin polarization of the tunneling electrons <cit.>.
It is determined by the spin dependence of the density of states near the Fermi energy of the ferromagnetic electrodes and the tunneling matrix elements for these electrons. 
Large TMR values at room temperature of up to several hundred percent<cit.> in magnetic-layer structured spin valves, also involving superconductors <cit.>, have been achieved. 
For low temperatures (T<5K), TMR values of about 1.9·10^4% have been observed in few-layer graphene heterotructures <cit.>. 
The TMR effect has already been utilized in magnetic random access memory (mRAM) devices which provide fast writing and reading as well as long lasting performance<cit.>.

Quantum dots placed between ferromagnetic leads allow to realize highly tunable spin-valve devices.
Control over single spins is possible by tuning gate voltages<cit.>.  
The Kondo effect and interaction-induced exchange fields have been measured<cit.>, and the compensation of the exchange field by external fields has been demonstrated <cit.>. 
Controlling the magnetoresistance through carbon-nanotube quantum dots by external gates experimentally has been reported <cit.>. 
Furthermore, spin precession in non-collinear quantum-dot spin valves has been tuned by electric gates <cit.>.

The theoretical description of the interplay between Coulomb interaction and spin-dependent tunneling in nonequilibrium situations imposed by a bias voltage is challenging. 
An analytical diagonalization of the underlying many-body Hamiltonian is impossible and controlled physical approximations are needed to investigate the emerging physical phenomena appropriately. Several theoretical approaches have been presented in recent years in order to describe quantum transport out of equilibrium in nanostructures.  
For weak tunnel couplings between quantum dot and leads, a perturbation expansion in the tunnel coupling is possible for various regimes of transport and gate voltages. 
A calculation to first order in the tunnel-coupling strength is sufficient to reveal the existence of an exchange field induced by finite Coulomb interaction on the quantum dot<cit.> or to show that the ferromagnetic Anderson-Holstein model acquires an effective attractive Coulomb interaction<cit.>. 
It has been shown that the TMR through a quantum-dot spin valve to first- plus second-order in tunneling features a zero-bias anomaly <cit.> that could be explained by the influence of sequential and cotunneling processes on accumulation and relaxation of the quantum spin.
A systematic scan of the TMR values (and their deviation from Julliere's value) as a function of gate and bias voltage has been performed within the same approximation scheme <cit.>.

With increasing tunnel-coupling strengths, however, higher-order tunneling processes become more and more relevant and perturbation theory fails.
Several numerical methods have been developed to include all orders in tunneling for the Anderson model.
For instance, numerical renormalization-group approaches have been used to address the Kondo problem in the absence of a bias voltage <cit.>.
Picking up the idea of renormalization, functional (fRG)<cit.> as well as time-dependent numerical renormalization-group theories (TD-NRG)<cit.> have been applied to the Anderson-Holstein model at finite bias voltages. 
For the interacting resonant-level model, a density-matrix renormalization group (DMRG) study has elaborated the low-energy physics and noise properties<cit.>. 
The time-dependent current through a dot featuring ringing effects have been studied within quantum-Monte-Carlo schemes (QMC)<cit.>. 
Furthermore, the Kondo regime of the Anderson model has been approached by multilayer multiconfiguration time-dependent Hartree theory<cit.>.

Some of these methods have already been used to investigate spin-dependent transport through a quantum-dot spin valve. 
It has been shown that within the linear-response regime NRG approaches work well<cit.>. 
Using a DMRG approach, the local density of states as well as the TMR has been calculated<cit.>. 
The Kondo problem for a quantum-dot spin valve in the nonlinear response regime has been addressed by an equations-of-motion method<cit.>. 
In addition, the influence of the interaction-induced exchange field on dipolar and quadrupolar spin moments has been investigated <cit.>. 

In this paper, we study the TMR of a quantum-dot spin valve in the resonant-tunneling regime and in the presence of both Coulomb interaction and a finite transport voltage.
For this, we adopt the numerically-exact method of iterative summation of path integrals (ISPI) <cit.> and generalize it to include spin-dependent tunneling.
The ISPI scheme is formulated on the Keldysh contour and builds upon a systematic truncation of real-time correlations induced by the leads.
In this paper, we concentrate on the tunneling current in the stationary limit, although time-dependent observables and higher-order correlation functions are easily accessible within the presented approach.
The TMR depends quite sensitively on the relative importance of the various tunneling processes.
In the high-temperature regime, our results coincide with a perturbative first- or first- plus second-order calculation.
With decreasing temperature, resonant tunneling becomes more and more important.
For weak to intermediate Coulomb interaction, we reach numerical convergence of our ISPI data down to low temperatures at which perturbation theory severely fails. 
Our method, thus, provides a suitable tool to study spin-dependent, non-equilibrium quantum transport in an interacting nanostructure, covering regimes that are inaccessible by other methods.

The structure of the article is as follows.
In Sec. <ref> we introduce the Hamiltonian of the quantum-dot spin valve and express the Keldysh partition function within a path-integral formulation.
The Coulomb term is treated with the help of a Hubbard-Stratonovich (HS) transformation that introduces an auxiliary, Ising-like spin field.
The ISPI scheme to integrate out the HS spins is described in detail in Sec. <ref>, which also contains checks of convergence and limiting cases.
In Sec. <ref>, we present our results and discuss the dependence of the TMR on temperature (<ref>) as well as on bias and gate voltage (<ref>). 
Finally, we conclude and give an outlook in Sec. <ref>.




§ MODEL & METHOD





 §.§ Hamiltonian


The Hamiltonian of a quantum dot coupled to two ferromagnetic leads (see Fig. <ref>) consists of three parts (we set ħ=1 throughout the paper), 

    ℋ = ℋ_dot + ℋ_leads + ℋ_T.

The first term 

    ℋ_dot = ∑_σ E_0 d^†_σ d_σ + U d^†_↑ d^†_↓ d_↓ d_↑

describes a single, spin-degenerate orbital of energy E_0, subject to Coulomb interaction U when the orbital is doubly occupied.
Operators d^†_σ/d_σ create/annihilate an electron with spin σ on the quantum dot, where σ =/ denotes the projection of the spin of the quantum-dot electron onto the quantization axis, which we choose along the magnetization direction of the source lead.
The energy E_0 is tunable via a gate voltage. 
To achieve a more compact notation, we combine the two annihilation operators d_σ to 𝒟 = (d_↑, d_↓)^T and similarly for the creation operators. 
After rewriting the interaction term by using the operator identity 2 d^†_↑ d^†_↓ d_↓ d_↑ = d^†_↑ d_↑ + d^†_↓ d_↓ -  (d^†_↑ d_↑ - d^†_↓ d_↓)^2, absorbing the quadratic part of the interaction into the single-particle energy, and defining ϵ_0 = E_0 + U/2, we arrive at 

    ℋ_dot = 𝒟^†ϵ_0σ_0 𝒟 - U/2(𝒟^†σ_z 𝒟)^2   .

Here, σ_0 is the 2× 2 identity and σ_z the z-Pauli matrix in spin space.

The electrons in the ferromagnetic leads α = L,R are described by the noninteracting Hamiltonian

    ℋ_leads = ∑_α𝐤𝒞^†_α𝐤ℰ_α𝐤𝒞_α𝐤,

with ℰ_α𝐤 = ([ ϵ_𝐤+ -μ_α         0;         0 ϵ_𝐤- -μ_α ]), where ϵ_𝐤τ is the single-particle energy for spin projection τ=± along the majority/minority spin direction, and μ_α denotes the chemical potential.
Again, we introduced spinors 𝒞_α𝐤 = (c_α𝐤+, c_α𝐤-)^T and similarly for 𝒞_α𝐤^†.
We assume the  density of states ρ_ατ=∑_ kδ(ω-ϵ_α kτ) to be constant in energy but spin dependent.
The degree of spin polarization is then characterized by p_α=(ρ_α +-ρ_α -)/(ρ_α ++ρ_α -), where p_α=0 corresponds to a nonmagnetic lead and p_α=1 represents a halfmetallic electrode with majority spins only.

Tunneling between quantum dot and leads is described by the Hamiltonian

    ℋ_T =    ∑_α𝐤𝒞^†_α𝐤 Y_α𝒟   + H.c.,

where t_α is the tunnel matrix element, Y_L=t_L σ_0 as well as Y_R=t_Rσ_0 for the parallel and Y_R=t_Rσ_x for the antiparallel magnetization configuration, respectively.

In this work we consider a symmetric setup, μ_L=-μ_R=eV/2, p_L=p_R=p, t_L=t_R=t, and a collinear magnetic configuration of the two leads.
Asymmetric choices for μ_α, p_α, t_α, or noncollinear magnetic configurations could straightforwardly be included in our framework.



 §.§ Keldysh functional integral

In the following we formulate the nonequilibrium theory for the quantum-dot spin valve in terms of coherent-state functional-integrals on the Keldysh contour. The Keldysh partition functional reads

    Z[η]=Tr[e^i(S_dot+S_leads+S_T[η])]   .

For each part of the Hamiltonian the respective action is given by

    S_dot =    ∫_𝒞 d t ψ̅_0(t) (i ∂_t - ϵ_0)σ_0 ψ_0(t)+ U/2(ψ̅_0(t)σ_z ψ_0(t))^2 
    
    S_leads =    ∫_𝒞 d t ∑_α𝐤ψ̅_α k(t) (i ∂_t σ_0 - ℰ_α𝐤)ψ_α k(t) 
    
    S_T[η] =    ∫_𝒞 d t ∑_α𝐤[ ψ̅_α𝐤 (t)Y_α(t) e^i η_α(t)ψ_0(t) + H.c.].

The path integral is calculated in the basis of fermionic coherent states, with ψ(t),ψ̅(t) denoting the corresponding Grassmann eigenvalue when acting with the annihilation/creation operator (spinors), as used e.g. in Hamiltonians Eqs. (<ref>) and (<ref>), onto a coherent state for lead and dot at time t along the Keldysh contour 𝒞, respectively<cit.>. 
We remark that the trace includes the Grassmann fields for both the quantum dot and the leads.
Since we have introduced a source term (counting field) η_α(t) to the tunneling action, we can obtain the current at time t_m via the functional  derivative

    I_α(t_m)=.-i eδ/δη_α(t_m)ln Z[η]|_η≡ 0.

Note that Eq. (<ref>) coincides with the expression for the current expectation value at measurement time t_m calculated from I_α(t_m)= e ⟨d/dt∑_𝐤𝒞^†_α k𝒞_α𝐤⟩(t_m).

The central quantity of interest for this paper is the tunneling magnetoresistance (TMR) defined by

    TMR≡I^α_p - I^α_ap/I^α_p  .

It measures the difference of the currents I^α_(p/ap) for a parallel (p), n̂_L=n̂_R, and an antiparallel (ap), n̂_L=-n̂_R, orientation of the leads' magnetization directions.
To obtain a dimensionless number between 0 and 1, we use as a normalization the current I^α_p for the parallel setup.
We remark that in many studies reported in the literature, the TMR is normalized with I^α_ap, which nominally gives a larger value of the TMR (can be even larger than 1) for the same physical effect.

Since the action depends quadratically on the Grassmann fields for the lead electrons, we can (without any approximation) perform the partial trace over the leads' degrees of freedom to obtain 

    Z[η]=Tr[ e^i(S_dot+S_env[η])]   ,

The trace is now over the dot degrees of freedom only, and the environmental action

    S_env[η]=∫_𝒞  dt ∫_𝒞 dt' ψ̅_0 (t)Σ_T(t,t',η)ψ_0(t')

involves the tunneling self energy in the presence of the counting field
Σ_T(t,t',η)=∑_αΣ^α_T (t,t',η). 
In the absence of interaction in the quantum dot, U=0, also the integration over ψ̅_0(t), ψ_0(t)-fields can be performed in Eq. (<ref>).
The resulting generating functional

    Z_0[η]=[D^-1(t,t',η)]   ,

is expressed in terms of the Green's function of the (noninteracting) quantum dot,

    D(t,t',η)=[(i∂_t-ϵ_0)σ_0-Σ_T(t,t',η)]^-1,

which accounts for tunneling via the tunneling self energy Σ_T(t,t',η).

We treat the Coulomb interaction term in the quantum-dot Hamiltonian by means of a Hubbard-Stratonovich (HS) transformation<cit.> for each time t along the Keldysh contour 𝒞.
To enable a numerical implementation, we first discretize the Keldysh contour into N time slices of length δ_t on the upper and N time slices on the lower contour, respectively.
As a consequence, the Green's function D becomes a (4N× 4N) matrix (the number of time slices N is to be multiplied with 2 due to spin σ= ↑, ↓ and another factor of 2 to distinguish the upper from the lower Keldysh contour). The HS transformation relies on the identity

    exp{-νiδ_t U/2( ψ̅^ν_0σ_z ψ^ν_0)^2}
    =1/2∑_s=±1exp(-sλ_νψ̅^ν_0σ_z ψ^ν_0 )

for the contribution of the Coulomb interaction term to the action at a given time slice of length δ_t on the upper/lower Keldysh contour, ν= ±. A solution for the HS parameter λ_ν is (for 0 < U < π/δ_t) given by<cit.>

    λ_ν=cosh^-1[cos(δ_t U/2)-iνsin(δ_t U/2)].

The main idea of the HS transformation is to replace a term that is quartic in the Grassman fields by terms that are only quadratic in the Grassman fields but coupled to a newly-introduced, Ising-like degree of freedom, s=± 1.
As a result, the quantum-dot fields can be integrated out, and the resulting Keldysh generating functional 

    Z[η]=∑_{ s}[Δ^-1[η]-Σ_C( s)] ,

becomes a multiple sum of determinants.
The sum runs over all 2^2N configurations of possible Ising-spin values for each time slice on the upper and lower Keldysh contour, combined in a multispin vector s = (s^+_1, s^-_1, s^+_2, s^-_2, …, s^+_N, s^-_N) of the HS parameters.

The matrix Δ[η] for the discretized Green's function of the noninteracting spin valve depends on the counting field η but not on the HS parameters s, while the matrix Σ_C( s) for the discretized self-energy due to Coulomb interaction is independent of the counting field but contains the HS parameters. 
Both are 4N× 4N matrices.
To specify their matrix elements, we make use of multi-indices a = (l,ν,σ), where the Trotter index l=1,…, N labels the time slice, the Keldysh index ν=± distinguishes the upper from the lower Keldysh contour, and σ=± labels the spin.
While the charging self energy is a diagonal matrix with matrix elements

    ( Σ_C( s) )_aa' = δ_ll'δ_νν'δ_σσ'σλ_ν s_l^ν,
 
the discretized Green's function Δ[η] is diagonal in spin but nondiagonal in the Trotter and the Keldysh indices.
In the absence of the counting field, we can make use of the time translational invariance and express the matrix elements of the Green's function (making use of the multi-indices a=(l,ν,σ) and a'=(l',ν',σ')) via its frequency representation,


    (Δ^-1[0])_aa' = δ_σσ'∫dω/2π e^-iω (l-l')δ_t[ (ω- ϵ_0) νδ_νν' 
    -i/2∑_αΓ_ασ( F_α (ω) )_νν'].


Here Γ_ασ=2 π|t_α|^2 ρ(ϵ_Fασ) characterizes the hybridization of the dot with lead α, and
( F_α (ω) )_νν' are the matrix elements of the 2× 2 Keldysh matrix 

    F_α(ω) =([  2f_α(ω)-1   -2f_α(ω); -2f_α(ω)+2  2f_α(ω)-1 ]).

The Fermi function f_α(ω)=[exp(β(ω-μ_α))+1]^-1 describes the equilibrium occupation distribution of lead α.
For later convenience, we define Γ_α=(Γ_α++Γ_α-)/2 as the average tunnel coupling for up- and down-spin electrons.
Furthermore, we assume a symmetric coupling to the left and right lead and define Γ=Γ_L=Γ_R.

In the presence of the counting field, time translation invariance is broken and there is no simple frequency representation of Δ^-1[η] similar to Eq. (<ref>).
However, since we are only interested in calculating the current (and not higher-order correlation functions), we can expand Δ^-1[η] to linear order in η_α at a single Trotter index m (that corresponds to the time t_m at which the current is measured) and drop the counting field anywhere else.
The resulting expression for Δ^-1[η] is given by Eq. (<ref>) with Γ_ασ being replaced by Γ_ασ[1- i δ_t/2η_α ( νδ_ml-ν' δ_ml')].

Furthermore, we observe that multiplying Δ^-1[η]-Σ_C( s) in Eq. (<ref>) with any matrix of equal dimension that depends neither on the counting field η nor on the HS parameters s will only result in multiplying Z[η] with an overall factor that drops out upon calculating the current  Eq. (<ref>).
We make use of this freedom to increase numerical stability by multiplying Δ[0].
So, in practice, we replace the Keldysh generating functional as defined in Eq. (<ref>) by

    Z[η]=∑_{ s} G   ,

with 

    G = [Δ^-1[η]-Σ_C( s)] ·Δ[0]   ,

where Δ^-1[η] has been expanded up to linear order in η_α(t_m) as explained above.




§ ITERATIVE SUMMATION OF PATH INTEGRALS




A brute-force summation over the 2^2N terms in Eq. (<ref>) is a hopeless task for realistic values of N.
A more sophisticated and efficient way to perform the sum is needed. 
Such a method is provided by the recently-established scheme of iterative summation of path integrals (ISPI) <cit.>. 
It is a completely deterministic approach, based on the iterative evaluation of traces over determinants as they appear in Eq. (<ref>).
The scheme relies on the fact that lead-induced correlations decrease with increasing time.
This allows us to reduce the number of summations from 2^2N to (N/K) 2^2K, where K is the number of Trotter slices (of size δ_t) over which correlations are taken into account.
This is a huge reduction for typical scenarios in which N is of the order of a few hundreds while K can be chosen to be 6 or 7.
We remark that a similar formulation in terms of the reduced density matrix of the system has been presented in Ref. Segal_2010.

In Fig. <ref>, we demonstrate how the dot's Green's function decays with time by showing the real and the imaginary part of one exemplary matrix element of Δ[0] as a function of t-t'=(l-l')δ_t in the continuum limit, δ_t → 0.
Increasing the gate and/or bias voltage, ϵ_0 and eV, respectively, may lead to faster oscillations but does not alter the exponential decrease<cit.>. 
This motivates us to neglect all contributions from time differences t-t' larger than a chosen time scale t_K.
In the discretized version this amounts to setting all matrix elements of Δ[0] with |l-l'| δ_t > t_K to 0. Therefore, lead-induced correlations are only included over a span of K Trotter slices, where K = t_K/δ_t.

The truncation after K Trotter slices makes the matrix G block tridiagonal.
It can, then,  be written in the form

    G= (
    [      G^11      G^12         0         0         ⋯         0;      G^21      G^22      G^23         0         ⋯         0;         0      G^32      G^33      G^34         ⋯         0;         0         0      G^43      G^44         ⋯         0;         ⋮         ⋮         ⋮         ⋮         ⋱         ⋮;         0         0         0         0         ⋯ G^N_K N_K ]).

The total 4N× 4N matrix is decomposed into blocks G^nn' of size 4K × 4K.
The index n runs from 1 to N_K=N/K, and only blocks G^nn' with |n-n'|≤ 1 have finite entries.
Since Σ_C( s) is a diagonal matrix, the definition of G in Eq. (<ref>) yields that the blocks G^nn' in the n-th row depend only on the 2K HS spins s_n=(s^+_(n-1)K+1, s^-_(n-1)K+1, …, s^+_nK, s^-_nK) from the Trotter slices (n-1)K+1 to nK but not from the ones outside this range.

To evaluate G for a given set of HS spins, we iteratively make use of Schur's formula 

    ( [ A B; C D ]) =  A · ( D - C A^-1 B )   .

In the first step, we set A = G^11.
To calculate the determinant of 4(N-K) × 4(N-K) matrix D - C A^-1 B, we again use Schur's formula with the new A taken as the upper left 4K × 4K block.
This procedure is repeated until we arrive at the lower right corner.
The result is

    Z[η]= ∑_{ s} G^11·Ǧ^22·…·Ǧ^N_KN_K

where Ǧ is iteratively defined by 

    Ǧ^nn = G^nn-G^n,n-1[Ǧ^n-1,n-1]^-1G^n-1,n.

together with Ǧ^11 = G^11.

Due to this iterative definition, Ǧ^nn depends on all the blocks G^mm' of Eq. (<ref>) that are placed above and/or to the left of G^nn, i.e. m,m' ≤ n.
This is, however, inconsistent with our decision to neglect all correlations beyond t_K.
Since G was approximated as being block tridiagonal, we should consistently eliminate the influence of G^mm' on Ǧ^nn when m,m' deviates from n by more than 1.
To accomplish this, we approximate Ǧ^nn by

    G̃^nn = G^nn-G^n,n-1[G^n-1,n-1]^-1G^n-1,n.

with G̃^11 = G^11, i.e., in comparison to Eq. (<ref>) we have replaced [Ǧ^n-1,n-1]^-1 by [G^n-1,n-1]^-1 on the right hand side.
Due to this approximation, G̃^nn and, therefore, also its determinant G̃^nn depends only on the 4K HS spins s_n and s_n-1 from the Trotter slices (n-2)K+1,…, nK but not from earlier or later ones.
We can arrange the 2^4K different values of G̃^nn for each HS spin configuration in a 2^2K× 2^2K matrix, the transfer matrix 

    Λ_n,n-1 = ( G̃^nn[𝐬_n, 𝐬_n-1] )

where each row corresponds to one of the 2^2K configurations for 𝐬_n and each column to one of the 2^2K configurations for 𝐬_n-1.
The Keldysh generating functional can, then, be written as a product

    Z[η] = 𝐞^T ·Λ_N_K,N_K-1·…·Λ_3,2·Λ_2,1·Λ_1,0

of these transfer matrices, where Λ_n,n-1 is a 2^2K× 2^2K matrix for n≥ 2, Λ_1,0 is a 2^2K dimensional column vector, and 𝐞^T = (1, …, 1) a 2^2K dimensional row vector.
Each of the N_K matrix multiplications requires a summation over 2^2K HS spin configurations. 
In total, there are N_K 2^2K summations instead of 2^2N without the truncation scheme.
(A more precise estimate of the scaling behavior of the numerical effort should take into account the effort for building G̃^nn from Eq. (<ref>).)

To illustrate the numerical effort, we mention that our code runs for given parameters and K=6 typically one hour on a linux cluster with 16 nodes when the current is measured at Γ t_m=10. For finer discretizations we use a CrayXT6m machine where data points including e.g. K=7  are calculated within 2.5 hours on 240 nodes. 



 §.§ Convergence


The discrete HS transformation Eq. (<ref>) needs finite discretization time steps δ_t.
This introduces a finite Trotter error for the calculation of the path integral along the Keldysh contour<cit.>, which scales with
δ_t^2 for δ_t→ 0.
A second source of systematic error comes from truncating correlations beyond t_K.
To reduce both errors and estimate the remaining error bars, we perform the following three-step procedure.

(i) First, we choose a finite correlation time t_K.
For this t_K, the time step is varied and we perform the well-established extrapolation<cit.> to the continuum limit, δ_t→ 0. 
A typical outcome of this procedure is shown in panel (a) of Fig. (<ref>).
We show the current for 𝐧̂_L=𝐧̂_R for different correlation times Γ t_K=0.8,1.0,1.2,1.5 and fixed Coulomb interaction U=Γ as a function of the Trotter step size. 
In this example we have chosen the other parameters as eV=ϵ_0=Γ, p=0.5, k_BT=0.2Γ. 
ISPI data are given as symbols, and the solid lines display the linear regression, which is used for the extrapolation towards δ_t→ 0.
This step eliminates the Trotter error.

(ii) We repeat the extrapolation δ_t→ 0 for different finite values of the correlation time t_K.
To eliminate the respectively associated truncation error, we extrapolate to infinite correlation time with the help of a linear regression for 1/(Γ t_K)→ 0.
This is illustrated in panel (b) of Fig. (<ref>) for different values of the Coulomb interaction U.

(iii) Finally, the error bars of the ISPI data are estimated from the standard deviations of both subsequent linear regressions. 

Since we are interested in the stationary values of the current (and the TMR), we need to set the measurement time t_m large enough to overcome transient behavior.
We find that Γ t_m≥ 10 is sufficient.

For the noninteracting case, U=0, the errors are eliminated very efficiently with the described extrapolations, such that error bars are smaller than the symbol size in all curves presented for this case.
With increasing U, the error bars become larger.
We are able to obtain reliable results within reasonable computation time for intermediate values of the interaction up to U=2 Γ. 
Whenever the orbital of the dot is tuned far outside the transport window opened by the bias voltage (large |ϵ_0| or small eV), the ISPI data for the current and the TMR become noisy.



 §.§ Benchmark 1: noninteracting case


To convince ourselves of the reliability of our ISPI code, we perform two benchmark tests.
First, we consider a noninteracting spin valve, U=0.
In this case, the current and, thus, the TMR can be calculated analytically.
The formulae for the current in the parallel and antiparallel configuration are given in Appendix <ref>.



In Fig. <ref>(a), we show the current through a noninteracting (U=0) quantum-dot spin valve for the parallel (blue) and the antiparallel (orange) configuration as a function of the bias voltage.
ISPI data are represented by the symbols, with error bars of the order of the symbol size and solid lines display the analytical results.
The parameters are chosen as ϵ_0=0, p=0.8, k_BT=0.2Γ.
We find excellent agreement for the full I(V) characteristics. For other parameters we have checked that the analytical curves are reproduced by the ISPI data with high precision as well.

The inset of Fig. <ref>(a) shows the TMR.
There is a well-pronounced zero-bias peak, which is well-resolved by the ISPI data.

In Fig. <ref>(b), we depict the TMR as a function of the level position ϵ_0 for different lead polarizations p. 
Again, the ISPI data (symbols) are in perfect agreement with the analytical solution (solid  lines). 
There is a pronounced peak at ϵ_0=0, again well-resolved by the ISPI data.
For increasing p, the current for the antiparallel lead configuration decreases, which results in a larger TMR, see Eq.(<ref>). 
 
 


 §.§ Benchmark 2: perturbation theory

 


In the presence of interaction, a full analytical solution is not available anymore.
In the limit of high temperature, however, the tunnel coupling can be treated in perturbation theory.
To lowest order, only sequential-tunneling processes contribute.
The first-order theory for a quantum-dot spin valve formulated within a diagrammatic real-time technique<cit.> has been extended to a full first- plus second-order (referred to as sequential plus cotunneling) calculation<cit.>.

This enables us to perform a second benchmark test.
We compare the ISPI data with those obtained from first- or first- plus second-order perturbation theory as a function of temperature and check whether the ISPI data converge to the perturbative results at high temperatures.
The result is shown in Fig. <ref>, where we plot the relative deviation

    Δ_I = | I^seq/cot-I^ISPI/I^ISPI|

of the first-order (seq) or first- plus second-order (cot) current from the ISPI results as a function of temperature.
For the perturbative results, we use the scheme developed in Ref. Weymann_2005.
We cover both the parallel and the antiparallel configuration.
Other parameters are set to eV=3Γ, ϵ_0=U=Γ, p=0.5 and lines are guides to the eye only.

We find the ISPI data in accordance with the perturbative calculation for k_BT≫Γ.
While a first-order description is only sufficient for temperatures much larger than Γ, the first- plus second-order calculation remains reliable for temperatures that are somewhat larger than Γ.
But for k_BT≲Γ, perturbation theory clearly fails and cannot compete with ISPI.



§ RESULTS


We now turn to the discussion of the TMR through a quantum-dot spin valve for the experimentally relevant regime in which the energy scales for Coulomb interaction, tunnel coupling, and temperature are all of the same order of magnitude.
Neither the fully analytical U=0 formulae, nor a perturbative treatment of tunneling, nor any zero-temperature calculation done with other methods such as NRG is applicable.

The TMR through a quantum-dot spin valve depends quite sensitively on the relative importance of different tunnel processes that contribute directly and indirectly to the current.
A reference is set by Julliere's value

    TMR^Jull = 2 p^2/1+p^2

describing a direct tunnel junction between two ferromagnets without a quantum dot placed in between.
(We remind that, in our definition (<ref>) of the TMR, we normalize with the current for the parallel configuration.
The more often used normalization with the current for the antiparallel configuration would give TMR^Jull = 2 p^2/(1-p^2).)
This value is easily derived from the fact that the transmission through the junction for a given spin is proportional to the corresponding density of states in the source and the drain, respectively, and the latter are proportional to 1+p for the majority and 1-p for the minority spins.
For the quantum-dot spin valve, the TMR is strongly influenced by the possibility to accumulate a finite spin polarization on the quantum dot.
The spin accumulation has the tendency to weaken the spin-valve effect, i.e., to reduce the TMR.
The value of the dot's spin polarization, on the other hand, results from an interplay of various processes that either lead to a spin accumulation or that provide a relaxation channel for the accumulated spin.
This makes the TMR a highly nontrivial function of temperature, gate- and bias voltage, and Coulomb interaction.
Since tunnel processes of different orders in the tunnel coupling contribute very differently to spin accumulation and spin relaxation, we are not surprised that a perturbative calculation of the TMR breaks down already at relatively high temperatures.

The fact that the TMR depends rather sensitively on the nature of the dominating transport channel has been extensively discussed in the literature.
Measurements of the TMR for weak tunnel coupling outside the Coulomb-blockade regime <cit.> could be explained within a sequential-tunneling picture.
In the Coulomb-blockade regime, however, cotunneling dominates, leading to an enhancement of the TMR as compared to the sequential-tunneling result <cit.>.
For a qualitative explanation of the TMR measured through a quantum-dot spin valve with strong tunnel coupling, a Breit-Wigner form of the transmission has been assumed <cit.>.
While the Breit-Wigner form takes resonant-tunneling processes into account, it fails to describe correlations that, at very low temperature, give rise to the Kondo effect.
The formation of Kondo correlations is accompanied with an enhancement of the transmission through the dot, which gives rise to a large, bias-voltage-dependent TMR as observed in Refs. Hauptmann_2008,Pasupathy_2004,Hamaya_2007.
Finally, we mention that also for quantum-dot spin valves in which not only the leads but also the island is ferromagnetic, an enhancement of the TMR due to cotunneling has been experimentally <cit.> and theoretically <cit.> found.




 §.§ Temperature dependence of TMR




The relative importance of resonant tunneling, i.e. higher-order tunneling contributions, as compared to sequential and cotunneling is nicely demonstrated by the temperature dependence of the TMR.
In Fig. <ref>, we show the results for a noninteracting quantum-dot spin valve, U=0, both (a) for the linear-response, eV=0.1 Γ, and (b) the non-linear response regime, eV= 3Γ.
The ISPI data (symbols) perfectly agree with the analytical results (ANA) available for the noninteracting case.
To distinguish the contributions from different orders in tunneling, we expand the analytical result up to order Γ, Γ^2, Γ^3, and Γ^4, respectively.

To lowest order, corresponding to the sequential-tunneling approximation, the TMR is temperature independent<cit.>. 
Including the second- (cotunneling), third- or fourth-order contribution makes the TMR a nonmonotonic function of the temperature.
Unphysical oscillations arise which, for small bias voltage, even extend to TMR values outside the interval [0;1].
The cotunneling approximation seems to give reliable results for temperatures somewhat smaller than k_BT≈ 10Γ.
Strikingly, even a fourth-order expansion of the full analytical result is severely failing to describe the correct temperature dependence of the TMR down to k_BT ≲Γ.
This demonstrates the overwhelming importance of resonant tunneling when temperature and Γ are of the same order of magnitude.
Only for k_BT≫ 10 Γ (not shown in the figure), the ISPI data converge to the sequential-tunneling result.



We proceed with presenting the temperature dependence of the TMR for an interacting quantum-dot spin valve in Fig. <ref>.
A full analytical solution to compare with is not available anymore.
We also refrain from showing results from a first- plus second-order calculation since the resulting TMR values are far off for the considered temperature regime.
Instead, we compare the TMR for different values of U.
We use the same parameters as for the non-interacting case and, again, show both the linear- and nonlinear-response regime.
The TMR increases with increasing Coulomb interaction U.

In the linear-response regime, Fig. <ref>(a), we find a non-monotonic dependence of the TMR on temperature.
For T→∞, the TMR approaches the sequential-tunneling value TMR=p^2=0.25, as required. 
Interestingly, for ϵ_0=0 and small bias voltage, also the low-temperature limit, T→ 0, yields TMR=p^2. 
In between, around k_BT ≈ 0.5 Γ, the TMR displays a local minimum with TMR ≈ 0.195. 

In Fig <ref>(b), we present the temperature-dependent TMR for the nonlinear bias regime, eV=3Γ.
In this case, the TMR increases monotonically with temperature and approaches the expected sequential-tunneling value p^2.



 §.§ Gate- and bias-voltage dependence of TMR




We proceed with analyzing the gate- and bias-voltage dependence of the TMR since both gate- and bias voltages are usually tuneable in experiment with high precision.
First, we consider the noninteracting case, U=0.
At low temperature, k_BT=0.1Γ, the TMR shows a rather rich structure as a function of both voltages, see Fig. <ref>.
In particular, there is a peak at ϵ_0=0, V=0.
The full width at half maximum of the peak as a function of ϵ_0 and V, respectively, depends on the degree of spin polarization p, see inset of Fig. <ref>.
With increasing bias voltage, the local maximum of the TMR as a function of ϵ_0 first stays at ϵ_0=0 but then splits into two local maxima.
In Fig. <ref>, the position of the local maximum as a function of ϵ_0 is highlighted by a dashed line.
Once the level position ϵ_0 is tuned far away from resonance, |ϵ_0|≫Γ, the TMR approaches Julliere's value 2p^2/(1+p^2)=0.4.
This is expected since in this regime, accumulation and relaxation of the quantum-dot spin can be neglected. 



To demonstrate once more the failure of perturbation theory, we show in Fig. <ref> the gate-voltage dependence of the TMR for U=0 in the linear-response regime, eV=0.1Γ, and compare with the full analytical result (ANA) as well as a perturbative expansion of the latter up to fourth order in Γ.
The ISPI data (symbols) perfectly match with the full analytic result, with error bars smaller than symbol size, while finite-order perturbation theory is not only quantitatively but even qualitatively wrong.
This emphasizes again the importance of resonant-tunneling effects for the TMR.



We now include finite Coulomb interaction U ≠ 0.
In Fig. <ref>(a), we show the gate-voltage dependence of the TMR for U/Γ=0, 1, and 2 in the linear-response regime, eV=0.1Γ. 
The local maximum at resonance, ϵ_0=0, remains clearly visible for all chosen values of U.
The error bars of the ISPI data are quite moderate even for U=2 Γ.
Since for large |ϵ_0|, the TMR goes up to Julliere's value, there are two minima symmetrically placed around the central maximum.
With increasing U, the depth of the minima seems to shrink a bit.

In panel (b) of Fig. <ref>, we switch to the nonlinear-response regime, eV=1.5Γ. 
This bias voltage is chosen close to the bifurcation point indicated in Fig. <ref>.
The central maximum has just split into two local maxima, a finite Coulomb interaction U seems to wash out the fine structure visible for U=0. 
Note that the limiting value of the TMR for |ϵ_0|≫Γ is again set by Julliere's value. 

We complete the discussion of the gate-voltage dependence of the TMR with the remark that also for the case of finite U, a first- plus second-order calculation completely fails at low temperatures.
This is shown in panel (c) of Fig. <ref>, in which the same parameters as in panel (a) are used.
 


Finally, we study the bias-voltage dependence of the TMR at ϵ_0=0 and k_BT=0.2Γ for different values of the Coulomb interaction strength, see Fig. <ref>.
Since the current is an antisymmetric function of V, the TMR is symmetric with respect to V=0. 
For large bias voltage, the TMR approaches the value obtained for sequential tunneling, TMR=p^2=0.25.

At U=0, there is a well-pronounced local maximum in the center.
Two local minima placed symmetrically around the central peak.
As U increases, the form of the TMR becomes modified.
At the center (linear-response regime), the TMR first grows until U≈Γ and, then, drops again.
In the non-linear regime, |eV|>Γ, the interaction always increases the TMR.
As a consequence, the minima are getting washed out for U> Γ.
The error bars indicate the increasing numerical challenge for increasing U, especially in the linear-response regime.



§ CONCLUSION AND OUTLOOK



The tunnel magnetoresistance through a quantum-dot spin valve is a highly non-trivial function of temperature as well as gate and bias voltage.
For negligible Coulomb interaction, the TMR can be calculated analytically.
In the limit of weak tunnel coupling, a perturbative treatment of tunneling including sequential or sequential plus cotunneling is possible.
In a realistic experimental scenario, however, the energy scales for charging energy, tunnel strength, and temperature are all of the the same order of magnitude, and neither neglecting interactions nor a perturbative treatment of tunneling is justified.
In contrast, resonant-tunneling effects in the presence of finite Coulomb interaction are important.

To tackle the task of calculating the TMR for experimentally relevant transport regimes, we make use of an iterative summation of path integral (ISPI) scheme, which is based on the evaluation of nonequilibrium path integrals on the Keldysh contour in an iterative and deterministic manner. 
We have generalized this scheme to include spin-dependent tunneling.
ISPI naturally includes tunneling contributions from all orders in the tunnel-coupling strength, which seems to be extremely important for a proper determination of the TMR.
For intermediate values of the charging energy, up to the same order of magnitude as the tunnel-coupling strength, numerical convergence is achieved.

After benchmarking against known results for limiting cases, we have analyzed the temperature as well as the gate- and bias-voltage dependence of the TMR.
We find that for k_B T ≲Γ, perturbation theory severely fails since resonant tunneling becomes important.
The TMR displays a well-pronounced peak at ϵ_0=0 and V=0. 
We are able to clearly resolve this peak within our ISPI scheme and to study the influence of a finite Coulomb interaction.

In conclusion, the ISPI treatment of the TMR through a quantum-dot spin valve offers a new theoretical tool that covers a transport regime in which other methods fail.
This includes the regime in which all the system parameters are of the same order and an expansion in a small parameter is impossible, which is particularly relevant in realistic experimental setups.




§ ACKNOWLEDGEMENTS

We acknowledge financial funding of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under project WE 5733/1-2, KO 1987/5-2 and – Projektnummer 278162697 – SFB 1242. Furthermore, we were inspired by fruitful discussions with Alfred Hucht on block matrix operations.



§ CURRENT FORMULAE FOR THE NONINTERACTING SPIN VALVE


In the absence of interaction, U=0, the current through a single-level quantum dot can be calculated analytically for all temperatures and bias voltages<cit.>.
Including a finite spin polarization of the leads is straightforward.
The current formulae for a parallel (p) and antiparallel (ap) magnetization configuration of the leads are


    I_p   =   ∫_-∞^∞ dω2  Γ^2_↓Γ^2_↑+(Γ_↓^2+Γ_↑^2) (ϵ_0 -ω )^2/(Γ_↓^2+(ϵ_0 -ω )^2) (Γ_↑^2+(ϵ_0 -ω )^2)[f(ω -eV/2)-f(ω+eV/2)]
    
    
    I_ap   =   ∫_-∞^∞ dω8Γ_↓Γ_↑/(Γ_↓+Γ_↑)^2+4(ϵ_0 -ω )^2[f(ω -eV/2)-f(ω + eV/2)]   .

These integrals are evaluated with the help of Cauchy's residual theorem.
Collecting the residues at the poles of the Fermi functions (given by the Matsubara frequencies) and of the expression in front of the Fermi functions, we get

    I_p   =   ∑_σ=↑,↓∑_γ=± 1∑_γ'=± 1γ/2Γ_σ[i γ' Ψ(1/2+γ'βΓ_σ/2 π -iβ/2π(γe V/2-ϵ_0))-π f(γe V/2 + iΓ_σ+ϵ_0)]
    
    
    I_ap   =   ∑_γ=± 1∑_γ'=± 1γΓ (1-p^2)[iγ' Ψ(1/2+γ'βΓ/2π-iβ/2π(γe V/2-ϵ_0))-π f(γe V/2+iΓ+ϵ_0)]   ,

where Ψ(…) denotes the digamma function.


50

Julliere_1975 M. Julliere, Tunneling between ferromagnetic films, Phys. Lett. A 54, 225-226 (1975).

Meservey_1994 R. Meservey, and P. M. Tedrow, Spin-polarized electron tunneling, Phys. Rep. 238, 173 (1994).

Parkin_2004 S. S. P. Parkin, Ch. Kaiser, A. Panchula, P. M. Rice, B. Hughes, M. Samant, and S.-H. Yang, Giant tunnelling magnetoresistance at room temperature with MgO (100) tunnel barriers, Nat. Mater. 3, 862 (2004).

Lee_2007 Y. M. Lee, J. Hayakawa, S. Ikeda, F. Matsukura, and H. Ohno, Effect of electrode composition on the tunnel magnetoresistance of pseudo-spin-valve magnetic tunnel junction with a MgO tunnel barrier, Appl. Phys. Lett. 90, 212507 (2007).

Ikeda_2008 S. Ikeda, J. Hayakawa, Y. Ashizawa, Y. M. Lee, K. Miura, H. Hasegawa, M. Tsunoda, F. Matsukura, and H. Ohno, Tunnel magnetoresistance of 604% at 300K by suppression of Ta diffusion in CoFeB/MgO/CoFeB pseudo-spin-valves annealed at high temperature, Appl. Phys. Lett. 93, 082508 (2008).

Vavra_2017 O. Vávra, R. Soni, A. Petraru, N. Himmel, I. Vávra, J. Fabian, H. Kohlstedt, and Ch. Strunk, Coexistence of tunneling magnetoresistance and Josephson effects in SFIFS junctions, AIP Advances 7, 025008 (2017).

Wang_2018 Z. Wang, I. Gutiérrez-Lezama, N. Ubrig, M. Kroner, M. Gibertini, T. Taniguchi, K. Watanabe, A. Imamoğlu, E. Giannini, and A. F. Morpurgo, Very large tunneling magnetoresistance in layered magnetic semiconductor CrI_3, Nat. Comm. 9, 2516 (2018).

Song_2018 T. Song, X. Cai, M. W.-Y. Tu, X. Zhang, B. Huang, N. P. Wilson, K. L. Seyler, L. Zhu, T. Taniguchi, K. Watanabe, et al., Giant tunneling magnetoresistance in spin-filter van der Waals heterostructures, Science 360, 1214 (2018).

Wolf_2010 S. A. Wolf, J. Lu, M. R. Stan, E. Chen, and D. M. Treger, The Promise of Nanomagnetics and Spintronics for Future Logic and Universal Memory, Proceedings of the IEEE 98, 2155 (2010).

Crisan_2015 A. D. Crisan, S. Datta, J. J. Viennot, M. R. Delbecq, A. Cottet, and T. Kontos, Harnessing spin precession with dissipation, Nat. Comm. 7, 10451 (2015).

Hauptmann_2008 J. R. Hauptmann, J. Paaske, and P. E. Lindelof, Electric-field-controlled spin reversal in a quantum dot with ferromagnetic contacts, Nat. Phys. 4, 373 (2008)

Pasupathy_2004A. N. Pasupathy, R. C. Bialczak, J. Martinek, J. E. Grose, L. A. K. Donev, P. L. McEuen, and D. C. Ralph, The Kondo Effect in the Presence of Ferromagnetism, Science 306, 86 (2004). 

Hamaya_2007
K. Hamaya, M. Kitabatake, K. Shibata, M. Jung, M. Kawamura, K. Hirakawa,  T. Machida, T. Taniyama, S. Ishida, and Y. Arakawa, Kondo effect in a semiconductor quantum dot coupled to ferromagnetic electrodes, Appl. Phys. Lett. 91, 232105 (2007).

Cottet_2006 A. Cottet, T. Kontos, S. Sahoo, H. T. Man, M.-S. Choi, W. Belzig, C. Bruder, A. F. Morpurgo, and C. Schönenberger, Nanospintronics with carbon nanotubes, Semicond. Sci. Technol. 21, S78 (2006).

Sahoo_2005 S. Sahoo, T. Kontos, J. Furer, C. Hoffmann, M. Gräber, A. Cottet, and C. Schönenberger, Electric field control of spin transport, Nat. Phys. 1, 99 (2005).

Koenig_2003 J. König and J. Martinek, Interaction-Driven Spin Precession in Quantum-Dot Spin Valves, Phys. Rev. Lett 90, 166602 (2003).

Braun_2004 M. Braun, J. König, and J. Martinek, Theory of transport through quantum-dot spin valves in the weak-coupling regime, Phys. Rev. B 70, 195345 (2004).

 Weiss_2015 S. Weiss, J. Brüggemann, and M. Thorwart, Spin-vibronics in interacting nonmagnetic molecular nanojunctions, Phys. Rev. B 92, 045431 (2015).

Weymann_2005 I. Weymann, J. König, J. Martinek, J. Barnaś, and G. Schön, Tunnel magnetoresistance of quantum dots coupled to ferromagnetic leads in the sequential and cotunneling regimes, Phys. Rev. B 72, 115334 (2005).

Weymann_2005_2 I. Weymann, J. Barnaś, J. König, J. Martinek, and G. Schön, Zero-bias anomaly in cotunneling transport through quantum-dot spin valves, Phys. Rev. B 72, 113301 (2005).

Costi_1994 T. A. Costi, A. C. Hewson, and V. Zlatić, Transport Coefficients of the Anderson Model via the Numerical Renormalization Group, J. Phys. Condens. Matter 6, 2519  (1994).

Khedri_2018 A. Khedri, T. A. Costi, and V. Meden, Nonequilibrium thermoelectric transport through vibrating molecular quantum dots, Phys. Rev. B 98, 195138 (2018).

Jovchev_2015 A. Jovchev and F. Anders, Quantum transport through a molecular level: a scattering states numerical renormalisation group study, Phys. Scr. T 165, 014007 (2015). 

Schmitteckert_2008 E. Boulat, H. Saleur, and P. Schmitteckert, Twofold Advance in the Theoretical Understanding of Far-From-Equilibrium Properties of Interacting Nanostructures, Phys. Rev. Lett. 101, 140601 (2008). 

QMC_RMP E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Continuous-time Monte Carlo methods for quantum impurity models, Rev. Mod. Phys. 83, 349 (2011).

Muehlbacher_Komnik_2008 T. L. Schmidt, P. Werner, L. Mühlbacher, and A. Komnik,  Transient dynamics of the Anderson impurity model out of equilibrium, Phys. Rev. B 78, 235110 (2008).

Thoss_2018 H. Wang and M. Thoss, A multilayer multiconfiguration time-dependent Hartree study of the nonequilibrium Anderson impurity model at zero temperature, Chem. Phys. 509, 13 (2018)

Weymann_2011 I. Weymann, Finite-temperature spintronic transport through Kondo quantum dots: Numerical renormalization  group study, Phys. Rev. B 83, 113306 (2011).

Simon_2007 P. Simon, P. S. Cornaglia, D. Feinberg, and C. A. Balseiro, Kondo effect with noncollinear polarized leads: A numerical renormalization group analysis, Phys. Rev. B 75, 045310 (2007).

Sindel_2007 M. Sindel, L. Borda, J. Martinek, R. Bulla, J. König, G. Schön, S. Maekawa, and J. von Delft, Kondo quantum dot coupled to ferromagnetic leads: Numerical renormalization group study, Phys. Rev. B 76, 045321 (2007)

Gazza_2006 C. J. Gazza, M. E. Torio, and J. A. Riera, Quantum dot with ferromagnetic leads: A density-matrix renormalization group study, Phys. Rev. B 73, 193108 (2006).

Swirkowicz_2008 R. Świrkowicz and M. Wilczyński and J. Barnaś, The Kondo effect in quantum dots coupled to ferromagnetic leads with noncollinear magnetizations: effects due to electron-phonon coupling, Journal of Physics: Condensed Matter 20, 255219 (2008).

Hell_2015 M. Hell, B. Sothmann, M. Leijnse, M. R. Wegewijs, and J. König, Spin Resonance without Spin Splitting, Phys. Rev. B 91, 195404 (2015).

Misiorny_2013 M. Misiorny, M. Hell, and  M. R. Wegewijs, Spintronic magnetic anisotropy, Nat. Phys. 9, 801 (2013).

Weiss_2008 S. Weiss, J. Eckel, M. Thorwart, and R. Egger, Iterative real-time path integral approach to nonequilibrium quantum transport, Phys. Rev. B 77, 195316 (2008); Phys. Rev. B 79, 249901(E) (2009).

Weiss_2013 S. Weiss, R. Hützen, D. Becker, J. Eckel, R. Egger, and M. Thorwart, Iterative path integral summation for nonequilibrium quantum transport, Phys. Stat. Sol. B 250, 2298 (2013).

Becker_2012  D. Becker, S. Weiss, M. Thorwart, and D. Pfannkuche, Non-equilibrium quantum dynamics of the magnetic Anderson model, New J. Phys. 14,  073049 (2012).

Huetzen_2010 R. Hützen, S. Weiss, M. Thorwart, and R. Egger, Iterative summation of path integrals for nonequilibrium molecular quantum transport, Phys. Rev. B 85, 121408(R) (2012).

Kamenev_2005 A. Kamenev, Many-body theory of non-equilibrium systems, in Nanophysics: Coherence and Transport, Les Houches Session LXXXI, edited by H. Bouchiat, Y. Gefen, S. Gueron, G. Montambaux, and J. Dalibard (Elsevier, New York, 2005)

Negele_Orland J. W. Negele and H. Orland, Quantum Many-Particle Systems, (Addison-Wesley, Redwood City, CA, 1988).

Hubbard_1959 J. Hubbard, Calculation of Partition Functions, Phys. Rev. Lett. 3, 77 (1959).

Hirsch_1983 J. E. Hirsch, Discrete Hubbard-Stratonovich transformation for fermion lattice models, Phys. Rev. B 28, 4059 (1983).

Segal_2010 D. Segal, A. J. Millis, and D. R. Reichman, Numerically exact path-integral simulation of nonequilibrium quantum transport and dissipation, Phys. Rev. B 82, 205323 (2010).

Weiss_2004 S. Weiss, and R. Egger, Path-integral Monte Carlo simulations for few electron quantum dots with spin orbit coupling, Phys. Rev. B 72, 245301 (2005).

Weymann_08 I. Weymann and J. Barnaś, Shot noise and tunnel magnetoresistance in multilevel quantum dots: effects of cotunneling, Phys. Rev. B 77, 075305 (2008).

Martinek_02 J. Martinek, J. Barnaś, S. Maekawa, H. Schoeller, and G. Schön, Spin accumulation in ferromagnetic single-electron transistors in the cotunneling regime, Phys. Rev. B 66, 014402 (2002).

Hamaya_2007_2
K. Hamaya, M. Kitabatake, K. Shibata, M. Jung, M. Kawamura, K. Hirakawa,  T. Machida, T. Taniyama, S. Ishida, and Y. Arakawa, Electric-field control of tunneling magnetoresistance effect in a Ni/InAs/Ni quantum-dot spin valve, Appl. Phys. Lett. 91, 022107 (2007).

Ono_97
K. Ono, H. Shimada, and Y. Ootuka, Enhanced Magnetic Valve Effect and Magneto-Coulomb Oscillations in Ferromagnetic Single Electron Transistor, J. Phys. Soc. Jpn. 66, 1261 (1997).

Schelp_97
L. F. Schelp,  A. Fert, F. Fettar, P. Holody, S. F. Lee, J. L. Maurice, F. Petroff, and A. Vaures, Spin-dependent tunneling with Coulomb blockade, Phys. Rev. B 56, R5747 (1997).

Yakushiji_01
K. Yakushiji, S. Mitani, K. Takanashi, S. Takahashi, S. Maekawa, H. Imamura, and H. Fujimori, Enhanced tunnel magnetoresistance in granular nanobridges, Appl. Phys. Lett. 78, 515 (1997).

Takahashi_98
S. Takahashi and S. Maekawa, Effect of Coulomb blockade on magnetoresistance in ferromagnetic tunnel junctions, Phys. Rev. Lett. 80, 1758 (1998).

Jauho_1994 A. P. Jauho, N. S. Wingreen, and Y. Meir, Time-dependent transport in interacting and noninteracting resonant-tunneling systems, Phys. Rev. B 50, 5528 (1994).





