Nonequilibrium Kondo effect in a magnetic field: Auxiliary master equation approach

We study the single-impurity Anderson model out of equilibrium under the influence of a bias voltage $\phi$ and a magnetic field $B$. We investigate the interplay between the shift ($\omega_B$) of the Kondo peak in the spin-resolved density of states (DOS) and the one ($\phi_B$) of the conductance anomaly. In agreement with experiments and previous theoretical calculations we find that, while the latter displays a rather linear behavior with an almost constant slope as a function of $B$ down to the Kondo scale, the DOS shift first features a slower increase reaching the same behavior as $\phi_B$ only for $|g| \mu_B B \gg k_B T_K$. Our auxiliary master equation approach yields highly accurate nonequilibrium results for the DOS and for the conductance all the way from within the Kondo up to the charge fluctuation regime, showing excellent agreement with a recently introduced scheme based on a combination of numerical renormalization group with time-dependent density matrix renormalization group.


Introduction
Since its discovery almost one century ago, the Kondo effect has been measured in many physical systems ranging from bulk materials to nanostructures. The latter are especially attractive to study, because the parameters controlling the effect can be precisely tuned in the laboratory. There is a variety of experiments on nanowires [1][2][3], two-dimensional electron gases confined in heterostructures [4,5], carbon nanotubes [6] and also organic molecules [7], to mention a few. Whereas a finite temperature and a bias voltage to probe the effect are perturbations that naturally arise in these experiments and should therefore be studied, it is also interesting to study the effect of an additional magnetic field.
It is known from these experiments that upon introducing a Zeeman magnetic field B the zero-bias conductance anomaly (i.e. the peak of the conductance G as a function of bias voltage f) splits into two peaks located at f increases almost linearly with B [1 -3, 8]. Theoretical calculations [2,[9][10][11][12][13][14][15]  , where T K is the Kondo temperature that characterizes the width of the zero-bias anomaly at zero temperature and zero field. At the same time, the magnetic field produces a similar split in the total impurity density of states (DOS) (spectral function), which again starts developing for magnetic fields of the order of the Kondo scale, and which corresponds to a shift B w  in the spin-resolved impurity DOS. However, in contrast to B f , this shift does not show the same strictly linear behavior. Accurate calculations based on Bethe ansatz and the numerical renormalization group (NRG) [16][17][18] [19]). Notice that less sophisticated equations of motion approaches [20] yield instead a constant slope of B w as well. On the other hand, the different behavior of B w and B f is in contradiction with the simple expectation [20] that the enhancement of the conductance should occur when the chemical potential difference reaches the splitting in the spectral function. Kondo physics out of equilibrium is a challenging issue from the theoretical point of view and it is hard to obtain accurate results for both the spectral function and the conductance for voltages beyond the linear-response regime, most nonequilibrium steady-state approaches being perturbative or their accuracy being uncontrolled.
In this paper, we investigate the single-impurity Anderson model (SIAM) in the presence of both a magnetic field B and a finite bias voltage f. We adopt the recently introduced auxiliary master equation approach (AMEA), which has been shown to produce very accurate results for spectral functions and current characteristics both in as well as out of equilibrium [21]. To confirm the accuracy of our results we compare them with the ones obtained within a hybrid method that combines NRG with the time-dependent density matrix renormalization group (tDMRG) [22] to address quantum impurities out of equilibrium. The two approaches compare excellently (see figure 6) also at zero bias voltage, where we directly compare the spectral function with NRG. Our results confirm the different behavior of B w and B f , showing that there is no incompatibility. We also evaluate the magnetization in the high and low field limit, confirming the presence of a plateau at high fields for bias voltages e g B B  f m | | observed in previous theoretical results [12]. This work is organized as follows: in section 2 the model and the solution method are described. We start with an introduction to the model, section 2.1, followed by a part about Keldysh Green's functions, section 2.2. Then the general idea of AMEA and the solution method are sketched, section 2.3. In section 2.4 the hybrid NRG-tDMRG method, which we use for comparison, is described. Section 3 contains the results and section 4 a summary and our conclusions.

Model
We study the SIAM in a magnetic field and out of equilibrium. Throughout this paper we use units of e k g 1

B B
 m = = = = | | and 1 G = , see equations (6) and (15) , U is the interaction strength and B the magentic field. The on-site energy f e s is chosen such that the system is particle-hole symmetric at B=0. The impurity is connected to two leads described by across the impurity. The Hamiltonian mediating the coupling between the impurity and the leads is given by with a symmetric hopping t t L R ¢ = ¢ . We assume that H leads produces a flat DOS r w l ( ) in the disconnected leads with a bandwidth of D 2 , where Θ is the Heaviside step function. In this flat-band model the hybridization strength Γ, defined in equation (15), is given by Using 1 G = as unit of energy yields t D ¢ = l p for the hopping to the leads. Throughout this paper we take D=10.
We furthermore use the following definition of the Kondo temperature T K , at B=0. G is the linear-response differential conductance, equation (18), ) .

Keldysh Green's functions
While there is only one independent Green's function in equilibrium, there are two in nonequilibrium: the retarded and the Keldysh Green's function, G R and G K , e.g., are independent of each other. At finite magnetic field they are furthermore different for both spin kinds. In steady state, when the system is time-translation invariant, they are defined as and in Fourier space, Upon introducing the Keldysh contour, these Green's functions can be arranged in a matrix structure, according to where the advanced Green's function is related to the retarded one by G G In this way, the familiar form of Dyson's equation is maintained, The hybridization function is given by where g w l ( ) is the (noninteracting) Green's function of the decoupled leads. Since these are in equilibrium, its components obey the fluctuation-dissipation theorem, ] denotes the Fermi function at temperature T and chemical potential m l . The DOS in the leads is connected to g R w l ( ), Therefore in equilibrium only one independent Green's function persists. The hybridization strength Γ is defined, using equation (12), Given the full interacting Green's function at the impurity, the spin-resolved and total spectral functions are calculated as The current across the impurity is determined via the Meir-Wingreen formula [20]. In case of a biasindependent lead DOS with L R r w r w = ( ) ( ), such as equation (5), it reduces to [23] j . In linear-response the differential conductance G j = f ¶ ¶ is calculated from equation (17) as is the Fermi function at zero bias. In the general case, we calculate the differential conductance from finite current differences using three-point Lagrange polynomials to approximate the derivative.

Method
We here present a short sketch of the AMEA used in this paper. For more details, we refer to [21,[24][25][26]. The idea is to map the physical system described by (1) to a finite and open auxiliary system that has almost the same hybridization at the impurity as the original one (12) and thereby maintains the impurity physics, which we are interested in. The auxiliary system consists of a small number of N B bath sites connected to Markovian environments and its dynamics is governed by a Lindblad master equation. The parameters in this equation are determined to achieve a corresponding auxiliary hybridization function aux w D ( )such that aux w w D » D ( ) ( )as accurately as possible, see [25]. The physical hybridization function D is calculated from the given lead DOS, equation (5), using equations (12)- (14) and the Kramers-Kronig relation that links the real and imaginary part of a Green's function. The auxiliary hybridization function aux D can be calculated for a general set of bath parameters by solving a noninteracting Lindblad problem, see, e.g. [24][25][26][27][28]. The determination of these parameters and thus the mapping to the physical system is carried out with a parallel tempering algorithm [25].
The resulting Lindblad equation is solved by using matrix product states (MPS) and the time evolving block decimation algorithm (TEBD), as described in [21]. Since the auxiliary Lindblad system is essentially exactly solvable, the approximation of the method lies in the difference between aux w D ( )and w D( ). As shown in [25], this difference vanishes exponentially upon increasing N B . Therefore, a moderate number of bath sites (N 14 20 B » -) is sufficient to reach the accuracy required in the present paper.
The results we present here are in the steady state, which is determined via time evolution and formally reached with t  ¥. The Green's functions are also calculated in the time domain, starting from the steady state; they are continued to large times by linear prediction and then subjected to a Fourier transformation. The bias voltage is realized by shifting the chemical potentials in the leads symmetrically with respect to each other,

Comparison to NRG-tDMRG quench calculations
We compare our data to results obtained in a hybrid NRG-tDMRG quench setup which is described in [22]. While AMEA treats the impurity model as a truly open quantum system in the sense of a Lindblad master equation, for 'small enough' time scales t one can equally well consider quenches in a closed quantum system [29,30]. Starting with an initial state in which the two leads are in thermal equilibrium, but held at different chemical potential, standard Hamiltonian time evolution will drive the system towards its 'steady state' until at some point in time finite-size effects set in. For the SIAM one faces the difficulty that the different energy and time scales inherent in the model have to be handled with care. The hybrid NRG-tDMRG approach presented in [22] meets this challenge by exploiting the fact that energy scales outside the transport window, where are effectively in equilibrium. Thus, they can be traced out using the NRG [31]. Subsequently, the non-equilibrium processes arising on the energy scale of the transport window are treated within this renormalized setup using a tDMRG [32][33][34][35] quench. Both methods, NRG and tDMRG, are implemented based on MPS.
For the high-energy range outside the transport window a logarithmic discretization is used, while the transport window itself is discretized linearly. After mapping the problem onto a chain, the Hamiltonian of the first part of this chain, which represents the high-energy modes, can be diagonalized using NRG. This yields a truncated effective low-energy basis for this part of the system, which can be seen as the local state space of a renormalized impurity (RI). This RI is coupled to the remainder of the leads, which corresponds to the energy range of the transport window and therefore has an effective bandwidth set by voltage and temperature. The quench is initialized with a state ini y Yñ = ñ Ä Wñ | | | , where ini y ñ | lies in the ground state sector of the RI and Wñ | is the thermal state of the remaining part of the leads at different chemical potential and decoupled from the RI. This state is time-evolved using tDMRG. The relevant time scale for this quench is given by the size of the transport window.
To further simplify the MPS calculation, the leads are described in the form suggested by the thermofield approach [36][37][38], in which the thermal state Wñ | is a pure quantum state, and, even more advantageously, a simple product state on the MPS chain. This implies that the time evolution of the tDMRG quench is started with a product state and, hence, with lowest possible entanglement.
In practice, the time evolution is typically limited, due to the entanglement growth, before finite size effects set in. So far, the approach has only been used to calculate expectation values, because the determination of spectral functions would need far more numerical resources. For all data points with D 0.14 f > there was no need to use NRG, because the transport window is of similar size as the full bandwidth. For high voltages convergence was achieved only in the current and not in the magnetization. However, the time dependence of the dot's occupation, n t á ñ s ( ), follows an exponential decay such that one can extrapolate to the steady-state value.

Results
Our approach allows for an accurate solution of the model in and out of equilibrium, below, but also above the energy scale T K , so as to take into account the influence of charge excitations and of the Hubbard bands. At the same time, below T K and in equilibrium our results show a remarkable agreement of the spectral function with NRG up to intermediate values of U 6  G , see figure 6(a) and [21]. Here we want to study the behavior and interplay of the spectral function and the differential conductance in the presence of a finite Zeeman magnetic field B and bias voltage f. In particular, we focus on the shift of the Kondo and of the zero-bias peak.
We start by plotting the impurity spectral function in equilibrium ( 0 f = ) for different magnetic fields B, see figure 1. Most of our results are obtained for an interaction of U=6, corresponding to a Kondo temperature of . At finite magnetic field, the spin degeneracy is lifted, resulting in different spectral functions for spin-up and spin-down electrons. At particle-hole symmetry they are related to each other, according to Upon increasing the magnetic field, the Kondo resonance is suppressed and it broadens, similarly to the effect of a bias voltage, see [14,21,24,[39][40][41][42][43][44][45][46][47][48][49]. Furthermore, a magnetic field causes a shift B w of the Kondo resonance to higher energies in the spin-resolved spectral function A  and produces a splitting A d in the total spectral function A A A ) . This splitting starts at B T K  , see also [11,17,50], and persists until the peaks merge with the Hubbard bands. The position of the Kondo resonance in A  becomes B » for large B, while for decreasing B the ratio B B w decreases (see figure 4), consistent with previous results, mainly on the Kondo model [9,10,[16][17][18]51]. Note that for large magnetic fields one has 2 , while for small magnetic fields A d is smaller, due to the overlap of the contributions from the two spin directions. A similar splitting is produced by a bias voltage in the absence of a magnetic field [14,21,24], so that it is interesting to study the combined behavior of the two effects. In the presence of both, a finite bias voltage and magnetic field, one would expect 4 peaks in the total spectral function at B 2 f   . This has been observed within an equation of motion approach in [10] (see also [12]). It is not easy to observe such a four-peak structure within a numerically controlled, nonperturbative approach. In our case, for U=6, the higher energy peaks merge with the Hubbard bands before the peaks are sufficiently far apart, so that they look more like shoulders than peaks. For this reason, we investigate this effect for U=8. Figure 2(a) shows the spin-resolved spectral functions A w  ( ) at B=1 for different bias voltages f and U=8. At 0 f = the position of the Kondo resonance B w is closer to B than for the U=6 case, due to the fact that T K is smaller here. As a result of the applied bias voltage the shifted Kondo resonance first acquires a broadening and then, starting from 1 f » , it gets split. The two peaks have a distance of f » as expected, but the splitting is not symmetric. The corresponding four-peak structure in the total spectral function can be seen in figure 2 [10]. A more direct quantity to be measured experimentally is the differential conductance G across the impurity. In figure 3 we plot the current j (a) as well as G (b) as a function of the bias voltage for different values of B. The parameters are the same as in figure 1. To test the approaches, in figure 6(b) we compare results from AMEA with the ones from the hybrid NRG-tDMRG calculation discussed in section 2.4. Results are essentially on top of each other. The magnetic field affects the zero-bias peak in the conductance by first broadening it up to B T K  and then producing a split [9-12, 14, 15], as observed experimentally [2,3,5,52]. Notice that G d , the splitting in G, starts at B 0.3 » and is slightly delayed in comparison to A d , the splitting in A w  ( ), figure 1(b), which sets in at B 0.2 » . The reason for the delay in the splitting is the averaging of the spectral function in the current integral (17), which smears out the effect of the split peaks. Since G G G = =   at particle-hole symmetry, B f , the shift in the spin-resolved conductance G  , exactly fulfills B 2 G f = d , in contrast to its spectral counterpart, On the other hand, the magnitude of the shift in G, while becoming B for B T K  , as shown in figure 4, it reaches this limit faster than the shift in A w  ( ). In fact, figure 4 suggests that, within the error bars 3 B f becomes B as soon as it shows up, in contrast to B w . This is consistent with experiments [1, 2, 52], which indicate a strictly linear behavior. At B f  but smaller than the bandwidth the differential conductance reaches a B-independent value of G G 0.27 0 » . Figure 5(a) shows the magnetization n n á -ñ   and 5(b) the double occupancy n n á ñ   at the impurity in dependence of the bias voltage for different magnetic fields. At large magnetic fields B T K  the magnetization shows a plateau for B  f followed by a logarithmic decrease (straight lines in figure 5(a)), in agreement with previous results, see [12]. At small magnetic fields B T K  it starts to decrease for T K f » . Again, we find a very good agreement between AMEA and NRG-tDMRG, see figure 6(c). For small magnetic fields the double occupancy has a minimum at 2 f » , which seems to be independent of T K , see [53]. This minimum vanishes at larger magnetic fields as the Zeeman splitting of the local level increases and hence presumably is governed by charge fluctuations. In figure 6 we display a comparison of results obtained within AMEA (dashed lines and circles) with results from NRG ((a) dotted lines) and the hybrid NRG-tDMRG scheme discussed in section 2.4 ((b), (c) squares). Equilibrium spectral functions (a), differential conductance (b) and magnetization (c) curves at different magnetic fields agree remarkably well between the two approaches. One can only see small deviations in the spectral functions at high energies, due to the logarithmic discretization in NRG, which makes it less accurate in this energy region. The inset in (a) shows a zoom around 0 w = , where NRG is known to produce essentially exact results. In this region the two spectral functions deviate by less than 1%. The differential conductance at finite bias, being evaluated from finite current differences (see section 2.2) in both approaches, is, in principle, more prone to errors. Nontheless, the results lie essentially on top of each other. On the other hand, as remarked in section 2.4, the magnetization from the NRG-tDMRG scheme is not fully converged to the steady state and the data have been extrapolated assuming an exponential decay of the occupancy n t á ñ s ( ). For this reason, at high voltages, we can see that the values for the magnetization lie slightly above the AMEA results. While it is, in principle, possible to calculate spectral functions within the NRG-tDMRG scheme, it is unclear at the moment, whether this is numerically feasible. For this reason, we do not provide a comparison between the two approaches in figure 6.

Summary and conclusions
In this paper, we studied the Anderson impurity model out of equilibrium under the influence of a bias voltage f and a magnetic field B. In particular, we addressed the issue of the different behavior of the shift of the Kondo   . Shift B f of the conductance peak (in figure 3(b)) and B w of the equilibrium spectral function (in figure 1(a)) divided by the magnetic field B plotted as a function of B. Parameters are as in figure 1. peak in the impurity spectral function and the one in the conductance anomaly as a function of the magnetic field. We also presented explicitly results for the spectral function showing a four-peak structure resulting from the combined effects of B and f.
Our results agree with previous theoretical and experimental results in the known limits B T K  and B T K  , while our approach allows us to access the intermediate regime B T , The key aspect of AMEA [21,[24][25][26] is that we can obtain very accurate results also for the spectral functions out of equilibrium, which is difficult by other methods. The accuracy of our results in the parameter regime we considered is confirmed by an excellent comparison of spectral functions with NRG at 0 f = (up to frequencies for which NRG is supposed to yield correct results), and of expectation values with a recently introduced hybrid NRG-tDMRG scheme [22] at finite bias voltages.
The two approaches adopted here, AMEA and the NRG-tDMRG scheme, deal with the challenge of describing the long time behavior of the nonequilibrium SIAM in a different manner. While AMEA explicitly describes an open quantum system and thus is not restricted to finite time scales, the quench approach renormalizes the problem down to the relevant energy scale. In addition, AMEA is able to evaluate the impurity spectral function. While, in principle, this is also possible in the NRG-tDMRG approach, from a numerical point of view, it would be more costly. Therefore, it is unclear at the moment, whether it is realizable in practice. Also for the magnetization AMEA was able to achieve better convergence, especially at high voltages.
In summary, it is convenient to use AMEA, whenever very long time scales are needed, or when information over the full energy range is required, as it is the case in the determination of spectral functions. For example, AMEA is an interesting tool for DMFT in nonequilibrium, where spectral functions are needed explicitly [26,[55][56][57][58][59]. On the other hand, the NRG-tDMRG approach is more flexible with respect to the parameter regime, as it uses an explicit renormalization of the impurity. In particular, it has proven to be able to describe very strong interactions such as U 12 G = and zero temperature T=0, see [22]. AMEA can deal with interactions of the same strength and temperatures down to T T 10 K [21]. Much larger values of U and/or much lower in T are not reachable at the moment, since we are limited in the number of bath sites 4 . This is also the reason, why we could not accurately check the well-known B T ln [19], in equilibrium. Our results may be consistent with a logarithmic asymptotics, but, in order to reliably confirm this behavior, we need to consider magnetic fields that are orders of magnitude larger and at the same time U  . Therefore, at the moment, it may be preferable to use the NRG-tDMRG quench approach, whenever it gets crucial to work in the scaling limit and for very low values of the bias voltage.
The only approximation in AMEA consists in replacing the physical bath hybridization function D with an auxiliary one aux D , so that the accuracy depends on the difference between the two functions. Of course, the corresponding error in the calculated results, e.g. the spectral function, is expected to be strongly frequency dependent, so that regions around the Fermi energies are probably more strongly affected. More specifically, due to the fact that at zero bias the Kondo scale depends exponentially on the 0 w = DOS, one may expect a corresponding exponential error in this scale. This is probably not yet the case at these moderate values of U 8  G used here, as can be deduced from our results in [21]. For larger U (and more bath sites), the way to avoid this exponential problem could be to carry out the fit by constraining