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

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–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 ϕ) splits into two peaks located at $\pm {\phi }_{B}$, where ${\phi }_{B}$ increases almost linearly with B [1–3, 8]. Theoretical calculations [2, 9–15] confirm this behavior showing an essentially constant slope, $e{\phi }_{B}\approx | g| {\mu }_{B}B$, almost all the way down to the point where the splitting disappears at $| g| {\mu }_{B}B\sim {k}_{B}{T}_{K}$, 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 $\pm {\hslash }{\omega }_{B}$ in the spin-resolved impurity DOS. However, in contrast to ${\phi }_{B}$, this shift does not show the same strictly linear behavior. Accurate calculations based on Bethe ansatz and the numerical renormalization group (NRG) [16–18] show that ${\omega }_{B}$ is initially smaller, starting as ${\hslash }{\omega }_{B}\approx \tfrac{2}{3}| g| {\mu }_{B}B$ and reaching $| g| {\mu }_{B}B$ for $| g| {\mu }_{B}B\gg {k}_{B}{T}_{K}$ (up to logarithmic corrections [19]). Notice that less sophisticated equations of motion approaches [20] yield instead a constant slope of ${\omega }_{B}$ as well. On the other hand, the different behavior of ${\omega }_{B}$ and ${\phi }_{B}$ 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 ϕ. 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 ${\omega }_{B}$ and ${\phi }_{B}$, 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\phi \lesssim | g| {\mu }_{B}B$ 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.

2.1. Model

We study the SIAM in a magnetic field and out of equilibrium. Throughout this paper we use units of ${\hslash }=e={k}_{B}={\mu }_{B}| g| =1$ and ${\rm{\Gamma }}=1$, see equations (6) and (15). The model is described by the following Hamiltonian,

$$\begin{eqnarray}&&H={H}_{\mathrm{imp}}+{H}_{\mathrm{leads}}+{H}_{\mathrm{coup}}.\end{eqnarray} \tag{ 1 }$$

${H}_{\mathrm{imp}}$ is the Hamiltonian of the impurity. It is a single-site Hubbard Hamiltonian with a spin-dependent on-site energy, accounting for the magnetic field,

$$\begin{eqnarray}&&{H}_{\mathrm{imp}}=\displaystyle \sum _{\sigma \in \{\uparrow ,\downarrow \}}{\varepsilon }_{f\sigma }^{}{f}_{\sigma }^{\dagger }{f}_{\sigma }^{}+{{Un}}_{f\uparrow }{n}_{f\downarrow },\end{eqnarray} \tag{ 2 }$$

with ${\varepsilon }_{f\sigma }=-\tfrac{1}{2}(U+\sigma B)$. ${f}_{\sigma }^{(\dagger )}$ is the fermionic annihilation (creation) operator at the impurity for spin σ, ${n}_{f\sigma }={f}_{\sigma }^{\dagger }{f}_{\sigma }^{}$, U is the interaction strength and B the magentic field. The on-site energy ${\varepsilon }_{f\sigma }$ is chosen such that the system is particle-hole symmetric at B = 0. The impurity is connected to two leads described by

$$\begin{eqnarray}&&{H}_{\mathrm{leads}}=\displaystyle \sum _{\lambda \in \{L,R\}}\displaystyle \sum _{k\sigma }{\varepsilon }_{\lambda k}^{}{d}_{\lambda k\sigma }^{\dagger }{d}_{\lambda k\sigma }^{}.\end{eqnarray} \tag{ 3 }$$

${d}_{\lambda k\sigma }^{(\dagger )}$ is the annihilation (creation) operator for electrons with spin σ in lead $\lambda \in \{L,R\}$ at level k (out of N energy levels); ${\varepsilon }_{\lambda k}$ is the energy of level k. The leads have different chemical potentials ${\mu }_{\lambda }$, realizing a bias voltage $\phi ={\mu }_{R}-{\mu }_{L}$ across the impurity. The Hamiltonian mediating the coupling between the impurity and the leads is given by

$$\begin{eqnarray}&&{{H}}_{{\rm{coup}}}=\displaystyle \frac{1}{\sqrt{{N}}}\displaystyle \sum _{{\lambda }\in \{{L},{R}\}}{{t}}_{{\lambda }}^{{\rm{{\prime} }}}\displaystyle \sum _{{k\sigma }}({{d}}_{{\lambda k\sigma }}^{\dagger }{{f}}_{{\sigma }}^{}+{\rm{H}}.{\rm{c}}.)\,\end{eqnarray} \tag{ 4 }$$

with a symmetric hopping ${t}_{L}^{{\prime} }{}_{}={t}_{R}^{{\prime} }{}_{}$. We assume that ${H}_{\mathrm{leads}}$ produces a flat DOS ${\rho }_{\lambda }(\omega )$ in the disconnected leads with a bandwidth of $2D$,

$$\begin{eqnarray}&&{\rho }_{\lambda }(\omega )=\displaystyle \frac{1}{2D}{\rm{\Theta }}(D-| \omega | ),\end{eqnarray} \tag{ 5 }$$

where Θ is the Heaviside step function. In this flat-band model the hybridization strength Γ, defined in equation (15), is given by

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{\pi }{2D}({{t}_{L}^{{\prime} }{}_{}}^{2}+{{t}_{R}^{{\prime} }{}_{}}^{2}).\end{eqnarray} \tag{ 6 }$$

Using ${\rm{\Gamma }}=1$ as unit of energy yields $t{{\prime} }_{\lambda }=\sqrt{\tfrac{D}{\pi }}$ 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,

$$\begin{eqnarray}&&G(T={T}_{K},\phi =0)=\displaystyle \frac{1}{2}{G}_{0},\end{eqnarray} \tag{ 7 }$$

at B = 0. G is the linear-response differential conductance, equation (18), ${G}_{0}=G(T=0,\phi =0)=1/\pi$.

2.2. 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

$$\begin{eqnarray}{G}_{\sigma }^{R}(t) & = & -{\rm{i}}{\rm{\Theta }}(t)\phantom{\rule{}{2.0ex}}\langle \phantom{\rule{}{2.0ex}}\{{f}_{\sigma }^{}(t),{f}_{\sigma }^{\dagger }\phantom{\rule{}{2.0ex}}\}\phantom{\rule{}{2.0ex}}\rangle ,\\ {G}_{\sigma }^{K}(t) & = & -{\rm{i}}\phantom{\rule{}{2.0ex}}\langle \phantom{\rule{}{2.0ex}}[{f}_{\sigma }^{}(t),{f}_{\sigma }^{\dagger }\phantom{\rule{}{2.0ex}}]\phantom{\rule{}{2.0ex}}\rangle ,\end{eqnarray} \tag{ 8 }$$

and in Fourier space,

$$\begin{eqnarray}&&{G}_{\sigma }^{\alpha }(\omega )=\int {G}_{\sigma }^{\alpha }(t)\exp ({\rm{i}}\omega t)\,{\rm{d}}{t},\end{eqnarray} \tag{ 9 }$$

with $\alpha \in \{R,K\}$. Upon introducing the Keldysh contour, these Green's functions can be arranged in a matrix structure, according to

$$\begin{eqnarray}{\underline{G}}_{\sigma }(\omega )=\left(\begin{array}{cc}{G}_{\sigma }^{R}(\omega ) & {G}_{\sigma }^{K}(\omega )\\ 0 & {G}_{\sigma }^{A}(\omega )\end{array}\right),\end{eqnarray} \tag{ 10 }$$

where the advanced Green's function is related to the retarded one by ${G}_{\sigma }^{A}(\omega )={G}_{\sigma }^{R}{(\omega )}^{\dagger }$. In this way, the familiar form of Dyson's equation is maintained,

$$\begin{eqnarray}{\underline{G}}_{\sigma }^{-1}(\omega ) & = & {\underline{g}}_{0\sigma }^{-1}(\omega )-\underline{{\rm{\Delta }}}(\omega )-\underline{{\rm{\Sigma }}}(\omega )\\ & = & {\underline{G}}_{0\sigma }^{-1}(\omega )-\underline{{\rm{\Sigma }}}(\omega ).\end{eqnarray} \tag{ 11 }$$

${\underline{G}}_{\sigma }(\omega )$ is the full interacting Green's function of the impurity connected to the leads, ${\underline{g}}_{0\sigma }(\omega )$ is the noninteracting Green's function of the disconnected impurity, $\underline{{\rm{\Delta }}}(\omega )$ is the hybridization of the impurity by the leads and $\underline{{\rm{\Sigma }}}(\omega )$ accounts for the interaction at the impurity. The noninteracting Green's functions are combined to ${\underline{G}}_{0\sigma }(\omega )={\underline{g}}_{0\sigma }^{-1}(\omega )-\underline{{\rm{\Delta }}}(\omega )$. The hybridization function is given by

$$\begin{eqnarray}&&\mathop{{\rm{\Delta }}}\limits_{\_}(\omega )=\displaystyle \sum _{\lambda }{{t}_{\lambda }^{{\prime} }{}_{}}^{2}{\mathop{g}\limits_{\_}}_{\lambda }(\omega ),\end{eqnarray} \tag{ 12 }$$

where ${\underline{g}}_{\lambda }(\omega )$ is the (noninteracting) Green's function of the decoupled leads. Since these are in equilibrium, its components obey the fluctuation-dissipation theorem,

$$\begin{eqnarray}&&{g}_{\lambda }^{K}(\omega )=2\pi {\rm{i}}(2{f}_{\lambda }(\omega ,T)-1){\rho }_{\lambda }(\omega ),\end{eqnarray} \tag{ 13 }$$

where ${f}_{\lambda }={[\exp [(\omega -{\mu }_{\lambda })/T]+1]}^{-1}$ denotes the Fermi function at temperature T and chemical potential ${\mu }_{\lambda }$. The DOS in the leads is connected to ${g}_{\lambda }^{R}(\omega )$,

$$\begin{eqnarray}&&{\rho }_{\lambda }(\omega )=-\displaystyle \frac{1}{\pi }{\mathfrak{I}}\,{g}_{\lambda }^{R}(\omega ).\end{eqnarray} \tag{ 14 }$$

Therefore in equilibrium only one independent Green's function persists. The hybridization strength Γ is defined, using equation (12),

$$\begin{eqnarray}&&{\rm{\Gamma }}=-{\mathfrak{I}}\,{{\rm{\Delta }}}^{R}(\omega =0).\end{eqnarray} \tag{ 15 }$$

Given the full interacting Green's function at the impurity, the spin-resolved and total spectral functions are calculated as

$$\begin{eqnarray}&&{A}_{\sigma }(\omega )=-\displaystyle \frac{1}{\pi }{\mathfrak{I}}\,{G}_{\sigma }^{R}(\omega ),\,\,\,\,\,{A}_{\uparrow \downarrow }=\displaystyle \frac{1}{2}({A}_{\uparrow }+{A}_{\downarrow }).\end{eqnarray} \tag{ 16 }$$

The current across the impurity is determined via the Meir–Wingreen formula [20]. In case of a bias-independent lead DOS with ${\rho }_{L}(\omega )={\rho }_{R}(\omega )$, such as equation (5), it reduces to [23]

$$\begin{eqnarray}&&j=\int {A}_{\uparrow \downarrow }(\omega )\,\gamma (\omega )({f}_{R}(\omega ,T)-{f}_{L}(\omega ,T)){\rm{d}}\omega ,\end{eqnarray} \tag{ 17 }$$

where $\gamma (\omega )=-{\mathfrak{I}}{{\rm{\Delta }}}^{R}(\omega )$. In linear-response the differential conductance $G=\tfrac{\partial j}{\partial \phi }$ is calculated from equation (17) as

$$\begin{eqnarray}&&G=\int {A}_{\uparrow \downarrow }(\omega )\,\gamma (\omega )\left(-\displaystyle \frac{\partial }{\partial \omega }f(\omega ,T)\right)\,{\rm{d}}\omega ,\end{eqnarray} \tag{ 18 }$$

where $f={f}_{L}{| }_{{\mu }_{L}=0}={f}_{R}{| }_{{\mu }_{R}=0}$ 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.

2.3. Method

2.4. Comparison to NRG-tDMRG quench calculations

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/{\rm{\Gamma }}\lesssim 6$, 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 ϕ. 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 ($\phi =0$) 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 ${T}_{K}\approx 0.2$. The temperature is fixed to $T/{T}_{K}\approx 0.25$. 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 ${A}_{\uparrow }(\omega )={A}_{\downarrow }(-\omega )$. 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–49]. Furthermore, a magnetic field causes a shift ${\omega }_{B}$ of the Kondo resonance to higher energies in the spin-resolved spectral function ${A}_{\downarrow }$ and produces a splitting ${\delta }_{A}$ in the total spectral function ${A}_{\uparrow \downarrow }=\tfrac{1}{2}({A}_{\uparrow }+{A}_{\downarrow })$. This splitting starts at $B\gtrsim {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}_{\downarrow }$ becomes $\approx B$ for large B, while for decreasing B the ratio ${\omega }_{B}/B$ decreases (see figure 4), consistent with previous results, mainly on the Kondo model [9, 10, 16–18, 51]. Note that for large magnetic fields one has ${\delta }_{A}=2{\omega }_{B}$, while for small magnetic fields ${\delta }_{A}$ is smaller, due to the overlap of the contributions from the two spin directions.

Zoom In Zoom Out Reset image size

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 $\pm B\pm \phi /2$. 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}_{\downarrow }(\omega )$ at B = 1 for different bias voltages ϕ and U = 8. At $\phi =0$ the position of the Kondo resonance ${\omega }_{B}$ 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 $\phi \approx 1$, it gets split. The two peaks have a distance of $\approx \phi$ as expected, but the splitting is not symmetric. The corresponding four-peak structure in the total spectral function can be seen in figure 2(b) with split peaks at $\omega \simeq \pm B\pm \tfrac{\phi }{2}$, see [10].

Zoom In Zoom Out Reset image size

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\gtrsim {T}_{K}$ and then producing a split [9–12, 14, 15], as observed experimentally [2, 3, 5, 52]. Notice that ${\delta }_{G}$, the splitting in G, starts at $B\approx 0.3$ and is slightly delayed in comparison to ${\delta }_{A}$, the splitting in ${A}_{\uparrow \downarrow }(\omega )$, figure 1(b), which sets in at $B\approx 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}_{\uparrow }={G}_{\downarrow }$ at particle-hole symmetry, ${\phi }_{B}$, the shift in the spin-resolved conductance ${G}_{\downarrow }$, exactly fulfills ${\phi }_{B}=\tfrac{{\delta }_{G}}{2}$, in contrast to its spectral counterpart, ${\omega }_{B}\geqslant \tfrac{{\delta }_{A}}{2}$. On the other hand, the magnitude of the shift in G, while becoming $\sim B$ for $B\gg {T}_{K}$, as shown in figure 4, it reaches this limit faster than the shift in ${A}_{\downarrow }(\omega )$. In fact, figure 4 suggests that, within the error bars³ ${\phi }_{B}$ becomes $\sim B$ as soon as it shows up, in contrast to ${\omega }_{B}$. This is consistent with experiments [1, 2, 52], which indicate a strictly linear behavior. At $\phi \gg B$ but smaller than the bandwidth the differential conductance reaches a B-independent value of $G\approx 0.27\,{G}_{0}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5(a) shows the magnetization $\langle {n}_{\uparrow }-{n}_{\downarrow }\rangle$ and 5(b) the double occupancy $\langle {n}_{\uparrow }{n}_{\downarrow }\rangle$ at the impurity in dependence of the bias voltage for different magnetic fields. At large magnetic fields $B\gg {T}_{K}$ the magnetization shows a plateau for $\phi \lesssim B$ followed by a logarithmic decrease (straight lines in figure 5(a)), in agreement with previous results, see [12]. At small magnetic fields $B\lesssim {T}_{K}$ it starts to decrease for $\phi \approx {T}_{K}$. 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 $\phi \approx 2$, 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 $\omega =0$, 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 $\langle {n}_{\sigma }\rangle (t)$. 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.

In this paper, we studied the Anderson impurity model out of equilibrium under the influence of a bias voltage ϕ and a magnetic field B. In particular, we addressed the issue of the different behavior of the shift of the Kondo 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 ϕ.

Our results agree with previous theoretical and experimental results in the known limits $B\ll {T}_{K}$ and $B\gg {T}_{K}$, while our approach allows us to access the intermediate regime $B,\phi \gtrsim {T}_{K}$ as well. The key aspect of AMEA [21, 24–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 $\phi =0$ (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–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/{\rm{\Gamma }}=12$ and zero temperature T = 0, see [22]. AMEA can deal with interactions of the same strength and temperatures down to $T\sim {T}_{K}/10$ [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⁴ . This is also the reason, why we could not accurately check the well-known $\sim {[\mathrm{ln}(B/{T}_{K})]}^{-2}$ behavior of ${A}_{\downarrow }(\omega =0)$ for $B\gg {T}_{K}$ [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 $\ll 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 $\underline{{\rm{\Delta }}}$ with an auxiliary one ${\underline{{\rm{\Delta }}}}_{{\rm{aux}}}$, 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 $\omega =0$ DOS, one may expect a corresponding exponential error in this scale. This is probably not yet the case at these moderate values of $U/{\rm{\Gamma }}\lesssim 8$ 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 ${\mathfrak{I}}{{\rm{\Delta }}}_{{\rm{aux}}}^{R}$ to coincide with ${\mathfrak{I}}{{\rm{\Delta }}}^{R}$ at $\omega ={\mu }_{R/L}$, or in any case require that the fit becomes more accurate around these points.

We would like to thank M Sorantin, I Titvinidze, W von der Linden and A Weichselbaum for fruitful discussions. DF, AD and EA were supported by the Austrian Science Fund (FWF) within the projects P26508 and F41 (SFB ViCoM), as well as NaWi Graz. FS and JvD were supported by the German-Israeli-Foundation through I-1259-303.10 and NIM. The calculations were partly performed on the D-Cluster Graz and on the VSC-3 cluster Vienna.