# Simulations of the dynamics of quantum impurity problems with matrix product
states

Matteo M. Wauters  Niels Bohr International Academy and Center for Quantum
Devices, Niels Bohr Institute, Copenhagen University, Universitetsparken 5,
2100 Copenhagen, Denmark  INO-CNR BEC Center and Dipartimento di Fisica,
Università di Trento, 38123 Povo, Italy  INFN-TIFPA, Trento Institute for
Fundamental Physics and Applications, Trento, Italy  Chia-Min Chung
Department of Physics, National Sun Yat-sen University, Kaohsiung 80424,
Taiwan  Center for Theoretical and Computational Physics, National Sun Yat-Sen
University, Kaohsiung 80424, Taiwan  Physics Division, National Center for
Theoretical Sciences, Taipei 10617, Taiwan  Lorenzo Maffi  Niels Bohr
International Academy and Center for Quantum Devices, Niels Bohr Institute,
Copenhagen University, Universitetsparken 5, 2100 Copenhagen, Denmark  Michele
Burrello  Niels Bohr International Academy and Center for Quantum Devices,
Niels Bohr Institute, Copenhagen University, Universitetsparken 5, 2100
Copenhagen, Denmark

###### Abstract

The Anderson impurity model is a paradigmatic example in the study of strongly
correlated quantum systems and describes an interacting quantum dot coupled to
electronic leads. Here we investigate its dynamics following a quantum quench
based on matrix product state simulations. We examine the behavior of its
impurity magnetization. Its relaxation allows us to extract the predicted
scaling of the Kondo temperature as a function of the impurity-lead
hybridization and quantum dot repulsion. Additionally, our simulations provide
estimates of the currents in the nonequilibrium quasi-steady state appearing
after the quench. Through their values, we examine the dependence of the
conductance on the voltage bias $V_{b}$ and on the impurity chemical potential
$V_{g}$, which displays a zero-bias Kondo peak. Our results are relevant for
transport measurements in Coulomb blockaded devices, and, in particular, in
quantum dots induced in nanowires.

##  I introduction

The Kondo effect is the most emblematic embodiment of strong correlations in
condensed matter systems. The advances in the fabrication and measurement
techniques of nanostructures allowed us to observe its distinctive zero-bias
conductance peak in a wide class of systems, including gate-defined quantum
dots [1, 2], nanotubes [3] and semiconducting nanowires [4].

In these mesoscopic systems, however, the dynamics of the quantum impurities
at the basis of the Kondo effect is typically too fast to be observed. A
complementary experimental platform has been recently offered by quantum
simulators of ultracold fermionic Yb atoms [5]. In these setups, the
characteristic time scales are much longer than in their solid state
counterpart, thus enabling the analysis of the dynamics of the spin impurities
at the basis of the Kondo effect in out-of-equilibrium transient states [6].

Inspired by these developments, in this work we analyze the dynamics of the
Anderson impurity model after a quantum quench through matrix product state
(MPS) simulations. By studying the transient behavior of its impurity
magnetization, we provide a numerical verification of the Kondo time scale
consistent with previous renormalization group results [7]. We derive the
conductance of the corresponding two-terminal problem by considering both
leads with a finite voltage bias and different charge impurity configurations.
Our results are relevant for the experimental study of Coulomb blockaded
nanowires with induced quantum dots.

The out-of-equilibrium properties of quantum impurity models following quantum
quenches are considered a paradigmatic playground to observe how strong
correlations develop through time evolution in many-body quantum systems and
have been recently studied by means of a vast set of analytical and numerical
techniques. The quench dynamics of the two-terminal Anderson impurity model
(AIM), in particular, has been addressed by quantum Monte Carlo methods to
estimate the related Keldysh Green’s functions [8, 9, 10, 11, 12, 13], time-
dependent numerical renormalization group [14], self-consistent diagrammatic
techniques [15], and tensor-network methods to study its real-time evolution
[16, 17, 18, 19, 20, 21, 22] (see also the comparison among different methods
in Ref. [23]).

Despite being one of the most studied interacting models in condensed matter
physics, the AIM keeps drawing attention as a standard test-bed for both
numerical and analytical investigations of out-of-equilibrium strongly
correlated systems and non-perturbative effects.

In the following we will benchmark the MPS algorithm we introduced in Ref.
[24] through its simulation of the time evolution of the two-terminal Anderson
impurity model (AIM). Our approach is tailored to examine the Kondo regime and
make predictions on its transport properties in the low-temperature limit.
Additionally, it can be readily extended to encompass more intricate forms of
hybridization between the leads and the impurity and characterize transport
features through more complex structures than the AIM, such as multilevel or
superconducting quantum dots and different kinds of interacting scatterers.

The rest of the paper is organized as follows. We introduce the model in Sec.
II and we give and overview of our numerical approach in Sec. III, summarizing
the main concepts of Ref. [24]. We illustrate our results in Sec. IV, where we
focus on the estimate of the Kondo temperature from the relaxation of the
magnetization and on the current-voltage curves. We summarize our work and
draw our conclusions in Sec. V. In the appendix A we make an explicit
comparison between the spin relaxation computed with our approach and the
results presented in Ref. [25] using the time-dependent numerical
renormalization group algorithm.

##  II The model

The AIM represents an electronic environment coupled with an interacting
magnetic impurity; it is one of the most popular yet simple models that
display the Kondo effect, and it constitutes the central element of dynamical
mean field theory methods for studying correlated materials, making it a
fundamental problem for many numerical algorithms [26].

Figure 1: (a) Sketch of the AIM: a single-level quantum dot with Hubbard
repulsion $U$ is tunnel-coupled to two noninteracting leads with chemical
potentials $\mu_{L}$ and $\mu_{R}$. (b) Schematic representation of the MPS
describing the system [24]. The sites of the chain represent single-particle
orbitals and are ordered by their energy. To account for the interaction, we
include an auxiliary bosonic charge site (represented by a square) which
counts the number of particles inside the dot. This construction introduces
long-range couplings (arrows) in the Hamiltonian MPO, which however do not
constitute an obstacle for the TDVP algorithm used for the time evolution.

Its Hamiltonian reads

$\widehat{H}=\widehat{H}_{\rm leads}+\widehat{H}_{\rm tunn}+\widehat{H}_{\rm
AIM}\ ,$ (1)

where

$\widehat{H}_{\rm
AIM}=U\hat{n}_{\uparrow}\hat{n}_{\downarrow}+V_{g}(\hat{n}_{\uparrow}+\hat{n}_{\downarrow})\,,$
(2)

with
$\hat{n}_{\sigma}=\hat{d}_{\sigma}^{\dagger}{\hat{d}^{\phantom{\dagger}}}_{\sigma}$
describing the occupation of the two spin states of a single-level quantum
dot, which, in turn, plays the role of the magnetic impurity and is
characterized by the Hubbard repulsion $U$ and the chemical potential $V_{g}$.
Unless otherwise stated, we will focus on the particle-hole symmetric point
$V_{g}=-0.5U$.

The lead Hamiltonian

$\widehat{H}_{\rm
leads}=-\sum_{\alpha,\sigma,l}t_{\alpha,\sigma,l}\left({\hat{c}^{\dagger}}_{\alpha,\sigma,l}{\hat{c}^{\phantom{\dagger}}}_{\alpha,\sigma,l}+{\rm
H.c.}\right)+\sum_{\alpha,\sigma,l}\mu_{\alpha,\sigma}\hat{n}_{\alpha,\sigma,l}$
(3)

describes two spinful fermionic chains ($\alpha=L,R$) with a spin-dependent
chemical potential $\mu_{\alpha,\sigma}$ and a hopping amplitude
$t_{l}=t_{0}\textrm{e}^{-(l-1)/\xi}$ that decreases exponentially as a
function of the distance from the site $l=1$, with a decay length $\xi$. This
is known as Wilson construction and it is commonly used in numerical
renormalization group approaches to impurity problems. Moreover, it has been
shown effective to increase both the resolution at small voltage bias, namely
by mimicking an effectively larger system, and the stability of the time
evolution in MPS simulation of transport problems [17, 24, 27]. Indeed, given
the finite size $L$, the density of states around the Fermi energy depends on
the hopping decay length $\xi$: the smaller $\xi$, the more states are shifted
toward the Fermi energy, leading to a smaller energy level spacing. Therefore,
a strong decay of the tunneling provides higher energy resolution to
accurately determine the dynamics for states close to zero energy (thus at
small bias voltages)[17].

Finally, the quantum dot and the leads are coupled with a standard tunneling
Hamiltonian

$\widehat{H}_{\rm
tunn}=-\sum_{\alpha,\sigma}J_{\alpha,\sigma}\left({\hat{c}^{\dagger}}_{\alpha,\sigma,1}{\hat{d}^{\phantom{\dagger}}}_{\sigma}+{\rm
H.c.}\right)\ ,$ (4)

where ${\hat{d}^{\phantom{\dagger}}}_{\sigma}$ destroys an electron with spin
$\sigma$ on the impurity level. Throughout this paper, we consider a uniform
tunneling strength between the quantum dot and the leads $J_{\alpha,\sigma}=J$
and we denote by $\Gamma=2J^{2}/t_{0}$ the effective tunneling rate in the
limit of infinite bandwidth (constant density of states).

To bring the system out of equilibrium, we adopt two different quantum quench
protocols [24, 28]: (i) in the zero-bias quench, we initialize the system with
$J=0$ and $\mu_{L}=\mu_{R}$, thus preparing a product state between the
impurity and the leads; at time ${\sf t}>0$, the leads are connected to the
quantum dot ($J>0$) and the system equilibrates towards a stationary state.
(ii) in the $\mu-$quench, the system is initialized in the ground state at
half filling, i.e. with uniform chemical potential $\mu_{L}=\mu_{R}$, and then
it evolves in time after a voltage bias $V_{b}$ is turned on. The first
protocol is more useful to study the relaxation of the impurity magnetization
and extract the Kondo temperature, while the second leads to a fast
convergence of the current to a nonequilibrium quasi-steady state (NEQSS) that
describes the dynamics of the quantum quenches of finite-size systems at
intermediate times. The time-evolution of the system after a $\mu-$quench is
indeed analogous to the dynamics of quantum quenches describing the melting of
domain walls in spin systems [29, 30] and is characterized by the onset of a
light-cone departing from the quantum dot at ${\sf t}=0$. The corresponding
post-quench Hamiltonian obeys a Lieb-Robinson bound [31, 32] such that, for
${t>0}$, the central region of the system (in proximity of the dot) is not
affected by finite-size effects until the quasi-particles produced by the
quench do not reach the boundary $l=L$ of the finite leads. As verified in
different models with fermionic leads, this implies that the conductance of
the system can be extracted by the currents in the NEQSS [33, 34, 24].
Furthermore, since in the $\mu$-quench the initial state already captures some
of the non-perturbative Kondo correlations, such state is closer to the Kondo-
like quasi-steady state that arises in transport measurements and it provides
the most efficient approach to capture the behavior of the currents in the
transient NEQSSs.

##  III Matrix product state simulations

Tensor networks offer a powerful framework to simulate the real-time evolution
of quantum impurity models [16, 17, 35, 18, 36, 37, 19, 38, 22, 20, 21, 39].
The MPS construction enhances the characterization of spectral features at
high frequencies, which hold significant importance in out-of-equilibrium and
time-dependent impurity problems [27].

In our simulations, we will describe the time-evolution of the closed and
finite system defined by the Hamiltonian (1) with leads of size $L$. As
mentioned above, we will focus our analysis to intermediate times after the
quantum quench, such that our results will be derived from NEQSSs which are
insensitive to finite-size effects. This implies, for instance, that,
differently from time-dependent NRG procedures, our MPS results are unaffected
by the incomplete thermalization caused by the discretized Wilson chain
discussed in Ref. [40] because they are unrelated to the thermal Gibbs state
emerging at infinite time.

Finite-size effects, however, set important limitations that we must consider.
First, in order to avoid boundary effects, our simulations must be restricted
to a finite time window set by the upper bound ${\sf t}_{\rm max}$ when the
quasiparticles produced during the quantum quench reach the edges of the
boundary and, for instance, cause a non-physical reflection of the current
after a $\mu$-quench. This problem is mitigated by the adopted Wilson chain
approach, which increases the typically accessible time-scale before observing
edge effects, ${\sf t}_{\rm
max}=\frac{\xi}{2t_{0}}\left(e^{\left(L-1\right)/\xi}-1\right)$ [24]. In the
uniform chain case instead, the current travels with constant speed
$v_{f}=2t_{0}$ through the leads, leading to an edge reflection after the time
interval $\left(L-1\right)/2t_{0}$. Additionally, finite-size effects
determine also an incomplete relaxation of the spin of the impurity (see next
section) after a zero-bias quench, which, again, sets a limit to the physical
time window we can consider. Finally, they also determine the resolution we
can achieve in the chemical potentials of the leads and, therefore, in the
introduction of a voltage bias in the system to estimate its conductance.

To simulate the post-quench dynamics, we model the system with the MPS
depicted in Fig. 1(b): each site represents a single-particle energy orbital
of the non-interacting and decoupled system ($U,J=0$), and we compute the
unitary time evolution of the closed system with the time-dependent
variational principle (TDVP) [41, 42].

In particular, we expand the construction presented in Ref. [24] with the
addition of the spin degrees of freedom; the MPS ”sites” are ordered based on
their energies [43], such that the entropy growth during the time evolution is
restricted in an energy window, thus a segment of the MPS, corresponding to
the voltage bias. This can be intuitively understood by considering that in
transport processes, only states in an energy window
$[\epsilon_{F}-eV_{b},\epsilon_{F}+eV_{b}]$, where $\epsilon_{F}$ is the Fermi
energy of the leads, contribute to the current. All the others remain in their
initial state, occupied or empty, and do not contribute to the entanglement
growth in the MPS. In particular, as pointed out in Ref. [43], the
entanglement entropy displays a logarithmic growth with time instead of a
linear increase. Since the basis states (MPS ”site”) are ordered by their
energies regardless the number of leads, introducing multiple leads (or the
spin degrees of freedom) is straightforward.

The interaction $U$ is introduced by including an auxiliary MPS site that
represents the charge $\widehat{N}=\hat{n}_{\uparrow}+\hat{n}_{\downarrow}$ of
the dot [24]. Tunneling events increase or decrease this charge by one. This
construction is not strictly necessary for a single impurity site as in the
AIM in Eq. (1), but it can easily allow for generalizations to multilevel dots
with a uniform all-to-all Coulomb repulsion described by an effective charging
energy. Moreover, representing separately the charge degree of freedom from
the dot particle number enables us to consider superconducting quantum
dots[24, 44].

In the chosen single-particle energy basis, $\widehat{H}_{\rm
leads}+\widehat{H}_{\rm AIM}(U=0)$ acts trivially on each site of the MPS as
they represent eigenstates of the Hamiltonian with $U,\ J=0$. Therefore the
dynamics is dictated only by the interplay between the tunneling Hamiltonian
coupling the leads with the quantum dot and the Hubbard repulsion
$U\hat{n}_{\uparrow}\hat{n}_{\downarrow}$. $\widehat{H}_{\rm tunn}$ is
nonlocal in this basis, but it can be described by a matrix product operator
(MPO) with limited bond dimension $\chi=8$, such that TDVP is not hampered by
the presence of these long-range tunneling elements and can be efficiently
used to simulate the dynamics for a long evolution time.

The method is implemented by using ITensor library [45]. The source code can
be found in Ref. 111The source code can be found in the repository
https://github.com/chiamin/QuenchTranportTwoChains.

##  IV Results

We first focus on the equilibration of the impurity after it is coupled to the
unbiased leads (zero-bias quench). In strongly correlated quantum impurity
models, the dynamics of the impurity magnetization is typically characterized
by two rates: $\Gamma$, which determines the short-time and nonuniversal
evolution; and the Kondo temperature $T_{K}$, whose inverse, the Kondo time
$\mathsf{t}_{K}=T_{K}^{-1}$, defines the time scale required for the formation
of the Kondo screening cloud (see, for instance, Ref. [25] for the Anderson
impurity model coupled with a single lead, and Ref. [47] for the corresponding
Kondo problem). In the renormalization group sense, the evolution for time
$\Gamma^{-1}<\mathsf{t}\lesssim\mathsf{t}_{K}$ is governed by the weak-
coupling fixed point of the Kondo problem [48], and $\mathsf{t}_{K}$
constitutes the decay time of the magnetization in this intermediate regime
towards the formation of a spin singlet with the conduction electrons.

Therefore, we aim to get an estimate of the Kondo temperature as a function of
the ratio between the interaction strength $U$ and the effective tunneling
rate $\Gamma$ from the dynamics of the impurity magnetization
$\langle\sigma^{z}\rangle=\langle\hat{n}_{\uparrow}\rangle-\langle\hat{n}_{\downarrow}\rangle$.
We prepare the quantum dot in the polarized state
$|\hat{n}_{{\uparrow}}=1,\hat{n}_{{\downarrow}}=0\rangle$ and measure its
evolution in time after a zero-bias quench. For this analysis, we choose
$L=64$ as the lead length and the hopping decay length between $\xi=8$ and
$\xi=32$, depending on the energy resolution needed to accurately measure the
magnetization up to times of the order of ${\sf t}_{K}$.

We consider two values for the interaction strength, $U=t_{0}$ and
$U=0.4t_{0}$, and we examine the particle-hole symmetric point $V_{g}=-0.5U$.
To extract the predicted exponential dependence of the Kondo temperature from
$U/\Gamma$ [49, 37, 48], we vary the hybridization strength $\Gamma$ between
$\sim U/20$ ($J\sim 0.15U$) and $\sim U/2$ ($J\sim 0.5U$).

Figure 2(a) shows the decay in time of the magnetization for different values
of $U/\Gamma$ while we fix $U=t_{0}$. We can easily identify three regimes: at
short time ${\sf t}\lesssim\Gamma^{-1}$, the different curves collapse on each
other as the relevant time scale for the relaxation of the impurity is set
only by $\Gamma$ (dashed black line). Indeed, notice that time is measured in
units of ${\Gamma^{-1}}$. At longer times, the relaxation rate depends on the
ratio $U/\Gamma$, with a slower decay the further the system lies in the
strongly interacting/weak-coupling regime. For these intermediate values of
$\mathsf{t}$, we can extract the relaxation time by exponential fits of the
data (dot-dashed gray lines). Finally, the impurity approaches a steady state
with a finite magnetization; in Fig. 2(a) this last regime is visible only for
$U=2\Gamma$. Due to the unitary dynamics, the system keeps the memory of its
initial state and a complete relaxation to a $SU(2)$ invariant state can not
be reached. Indeed, the initial state is prepared with unpolarized leads and
the impurity in the $|\uparrow\rangle$ state, leading to a total magnetization
of 1. as the system relaxes, we expect a residual impurity magnetization
$\langle\sigma^{z}\rangle\sim(2L)^{-1}$, assuming it spreads over the whole
system. Comparable results have been obtained for the AIM with a single lead
in Ref. [37] through the real-time density matrix renormalization group (DMRG)
applied to a superposition of four MPS which results in a logarithmic growth
of the entanglement similarly to our energy-basis construction.

Figure 2(b) illustrates the inverse of the relaxation times
$\mathsf{t}_{K}(U/\Gamma)$ extracted from the magnetization decay at
intermediate times [grey lines of panel (a)] as a function of $U/\Gamma$ and
for two values of the Hubbard interaction $U$. We interpret this quantity as
the Kondo temperature $T_{K}\sim\mathsf{t}_{K}^{-1}$. A comparison with the
renormalization group prediction

$T_{K}\sim\sqrt{U\Gamma}\textrm{e}^{-\frac{\pi U}{8\Gamma}}$ (5)

(solid black line) shows excellent agreement with our data for both values of
the interaction strength. In the inset, the same data are displayed in
logarithmic scale to emphasize the exponential dependence of the Kondo
temperature from $U/\Gamma$. Moreover, the two datasets perfectly collapse on
top of each other, highlighting the universal character of the exponential
decay linked to the Kondo temperature. At very weak coupling ($U/\Gamma\gg 1$)
the long evolution time needed for an accurate estimate of the relaxation time
can not be reached and our data deviate from the analytical prediction.

Figure 2: (a): Expectation value of the magnetization
$\langle\sigma^{z}\rangle=\langle\hat{n}_{\uparrow}\rangle-\langle\hat{n}_{\downarrow}\rangle$
as a function of time (in units of $\hbar/\Gamma$) in a zero-bias quench. The
dashed black lines indicates the short-time relaxation $\textrm{e}^{-{\sf
t}\Gamma/\hbar}$ while the gray dot-dashed lines highlight the slow dynamics
due to the Kondo resonance $\textrm{e}^{-{\sf t}/{\sf t}_{K}}$. The data
correspond to $U=t_{0}$ (b): Kondo temperature, extracted from the slow
relaxation shown in panel (a), versus the effective interaction strength for
two values of $U$. The solid black line corresponds to the RG prediction for
the Kondo temperature $T_{K}\sim\sqrt{U\Gamma}\textrm{e}^{-\frac{\pi
U}{8\Gamma}}$. The inset shows the same data in logarithmic scale to emphasize
the exponential behaviour of $T_{K}$. For large values of $U/\Gamma$ our
accuracy is limited by the small signal-to-noise ratio of the time evolution
of the magnetization.

Although the entanglement growth ultimately limits our ability to simulate the
evolution of large systems at long time, thus preventing the observation of
Kondo correlations for very weak coupling, our method allows us to observe
nonperturbative effects emerging directly from the nonequilibrium properties
of the AIM.

Figure 3: Current vs voltage bias in the symmetric point $V_{g}=-0.5U$ for
three different values of the hybridization strength $\Gamma$ and $U=t_{0}$.
Our data (markers) well match NRG results (solid lines); the temperature of
the NRG calculation is $T=10^{-6}U$. The dashed line corresponds to the
quantized current $I=2\frac{e^{2}}{h}V_{b}$.

Figure 4: (a) Differential conductance as a function of the induced charge
$n_{g}$ and of the voltage bias $V_{b}$ between the left and right leads, in
the strongly interacting/weak-coupling regime $U/\Gamma=12.5$. The zero bias
peak extending between the two sequential tunneling resonances at $n_{g}=0.5$
and $n_{g}=1.5$ signals the onset of the Kondo effect. (b) Differential
conductance at zero bias as a function of $n_{g}$, extrapolated with finite
bias current obtained at $V_{b}=0.04U$. The color legend is the same as of
Fig. 3.

The data in Fig. 2 are obtained with a zero-bias quench, i.e., with a
vanishing bias voltage. When dealing with transport properties, we evolve
instead the system with a voltage bias $V_{b}=\mu_{L}-\mu_{R}$ between the two
leads in order to observe a quasi-stationary current. Here we use the
$\mu-$quench protocol: in this scenario the initial state is the correlated
ground state of the Hamiltonian $\hat{H}$ (obtained through the DMRG),
quenched at ${\sf t}=0$ to a Hamiltonian with a finite voltage bias. To
simulate the quench dynamics at finite bias, we also need to adjust the decay
length of the hopping amplitude in the lead, $\xi$, such that the density of
states in the leads is approximately constant in the energy interval between
$\mu_{L}$ and $\mu_{R}$. The convergence in the simulation parameters ($\xi,\
L$, and the TDVP time step discretization) is reached when the current signal
displays a plateau in time long enough to reliably extract its expectation
value in the NEQSS that develops after the quantum quench. The maximum bond
dimension adopted is $\chi=2000$ with a truncation error $O(10^{-8})$.

In Fig. 3 we plot the quasi steady current as a function of the voltage bias
for different values of the effective tunneling rate $\Gamma$, while keeping
fixed the Hubbard interaction $U$, at the particle-hole symmetric point
$V_{g}=0.5U$. As we approach the strong-coupling regime $\Gamma\sim U$, the
current tends toward a linear response with a quantized differential
conductance $\frac{\mathrm{d}I}{\mathrm{d}V_{b}}=2\frac{e^{2}}{h}$, i.e. there
are two perfectly transmitting channels (dashed black line). The Kondo
temperature $T_{K}$ sets the extension of the bias window in which this
quantization occurs [50]. In particular, as showed in Fig. 2(b), the Kondo
energy scale drops down exponentially at weak coupling, $U/\Gamma\gg 1$, and
can become smaller than the values of voltage bias we can resolve with the
chosen lead length $L=100$ and hopping decay length $\xi=30$. This explains
the apparent deviation from the quantized conductance at weak coupling
$U/\Gamma=12.5$ in Fig 3. Our results well match the estimate of the current
using the Landauer-Büttiker formula and the spectral function of the AIM
obtained from NRG calculation [51, 52] (solid lines). Analogous results have
been recently obtained in Ref. [22] based on an MPS representation, in the
temporal domain, of the leads described as Feynman-Vernon influence
functionals. We remark that, away from the strong coupling regime, we can
simulate transport for voltages larger than the tunneling rate $\Gamma$. The
main limitation comes from the faster entanglement growth when states in a
large energy window contribute significantly to transport, which happens when
$V_{b}$ covers a significant fraction of the leads’ bandwidth.

Figure 4(a) illustrates the differential conductance in the weak-coupling
regime ($U/\Gamma=12.5$) as a function of the bias $V_{b}$ and the induced
charge parameter $n_{g}$ which is linked to the chemical potential as
$V_{g}=\frac{U}{2}(1-2n_{g})$ and determines the expectation value of the
total occupation of the quantum dot. We derive the differential conductance in
Fig. 4(a) from the simulation of a $\mu$-quench protocol in which the system
is initialized in the ground state of $\widehat{H}$ at half filling, thus for
the particle-hole symmetric point ($n_{g}=1,V_{g}=-0.5U$). At time
$\mathsf{t}=0$, both the induced charge $n_{g}$ and the bias voltage $V_{b}$
are quenched to their final value [horizontal and vertical axis of Fig. 4(a)].

At $n_{g}=0.5$ and $n_{g}=1.5$, we observe two bright zero-bias sequential
tunneling resonances, corresponding to the degeneracies between the empty and
singly-occupied dot $(n_{g}=0.5)$ or between the singly and doubly occupied
dot $(n_{g}=1.5)$. At finite voltage, the conductance peaks are prolonged
along the lines $V_{b}=\pm U(1-2n_{g})$ and $V_{b}=\pm U(3-2n_{g})$, following
the resonances between each biased lead and the quantum dot. Between the two
charge-degeneracy points, an extended zero bias peak indicates the onset of
the Kondo effect, although, for strong interaction and weak coupling, we can
not see the quantization of the conductance. This limitation originates mainly
from the high voltage resolution needed to sample the current at energies
below the Kondo temperature, which for $U/\Gamma=12.5$ is of the order
$T_{K}\sim 10{{}^{-3}}U$, where we expect the quantized linear response. To
reach such a resolution in $V_{b}$, we need either a larger system size or a
shorter decay length $\xi$. The former makes the simulations computationally
more expensive, while the latter induces nonphysical effects in simulations at
higher energy, preventing the calculation of the differential conductance in a
wide bias range. As common in nanostructure experiments (see, for instance,
Ref. [4]), this zero-bias peak does not extend to $n_{g}<0.5$ or $n_{g}>1.5$
where the ground state of the quantum dot becomes, respectively, empty or
fully occupied, thus losing the doublet degeneracy necessary for the Kondo
effect.

This last feature is particularly evident by looking at the zero-bias
conductance as a function of $n_{g}$, shown in Fig. 4(b). For all coupling
strengths, the conductance is suppressed in the absence of the doublet
degeneracy ($n_{g}<0.5$ or $n_{g}>1.5$) in the ground state. For
$U/\Gamma=3.13$ (green triangles) the quantized peak
$\frac{\mathrm{d}I}{\mathrm{d}V_{b}}=2\frac{e^{2}}{h}$ is reached, as
indicated by the dot-dashed grey line. As $\Gamma$ is decreased, the finite
bias used to compute the conductance becomes too large to correctly capture
the quantization in the linear response regime but the enhancement of the
conductance due to the Kondo resonance remains clearly visible for $0.5\leq
n_{g}\leq 1.5$. Figure 4(b) also exemplifies the expected crossover from a
weak-coupling to a strong-coupling regime of the linear conductance extracted
from the current at finite DC voltage, as the tuning of $\Gamma$ brings the
system deeper into the Kondo regime.

##  V Conclusions

In this work, we applied the tensor network method introduced in Ref. [24] to
study the Kondo effect in the Anderson impurity model. In particular, we used
an MPS+TDVP approach to study the dynamics of a single-level interacting
quantum dot coupled to two fermionic leads after quantum quenches of the
Hamiltonian parameters. We examined both the out-of-equilibrium evolution of
the quantum dot magnetization and the electric transport features emerging in
a nonequilibrium quasi-steady state after the quench.

The magnetization dynamic allows us to obtain a good estimate of the Kondo
temperature as the inverse of its relaxation time when the quantum dot is
coupled to unbiased leads. Such estimate is in agreement with renormalization
group results [7]. In particular, the magnetization decay displays two typical
time scales: the effective coupling rate with the leads and the Kondo time
scale. The appearance of these two decay regimes for short and intermediate
times is reminiscent of the experimental results concerning the evolution of
the spin population of impurities in 1D ultracold Yb gases [5].

Concerning the study of the conductance of the system, relevant for transport
measurements in nanostructures, our simulations allow us to study its
evolution when a voltage bias is applied between the two leads. By looking at
the emergent quasi-steady state, we can reconstruct its Coulomb blockade
properties as well as the emergence of a Kondo peak at zero bias. The latter
appears when the impurity chemical potential fixes its ground state in the
degenerate singly-occupied sector and the related differential conductance
approaches the quantized value $G=2\frac{e^{2}}{h}$ in the strong-coupling
regime.

We can simulate the system dynamics in a broad parameter range, from a
strongly interacting/weak-coupling regime to a strong-coupling one, well
beyond the applicability of standard perturbative master-equation approaches.
Moreover, our method is not limited by single-site or small impurities but can
be easily extended to multilevel quantum dots or nanowires with long-range
Coulomb repulsion. Despite using Wilson chains to model the leads, our
approach does not suffer from problems due to their finite heat capacity [40].
Indeed, our MPS simulations extract transport properties from intermediate-
time dynamics where the energy introduced after the quench does not yet
translate into a sizeable effective temperature, differently from the eventual
asymptotic steady state.

Additionally, our approach can address superconducting systems, opening the
path for the study of the out-of-equilibrium dynamics of the topological Kondo
effect [53, 54, 55], which arises in multiterminal impurities with p-wave
superconducting coupling (see Ref. [44]). This kind of system can be easily
described by identifying the spin degrees of freedom of the AIM as a label for
different spinless leads.

In general, our method can thus be used to examine transport phenomena in
hybrid superconducting-semiconducting multiterminal devices with strong
Coulomb interactions (see, for instance, Refs. [56, 57]), without being
limited to a weak-coupling regime. This offers the possibility of
investigating the variety of subgap states [58] that can appear in these
platforms, thus providing important details towards the realization of
Majorana - Cooper pair boxes and other building blocks for quantum devices
[59, 60].

## Acknowledgements

We thank J. Paaske and V. Baran for fruitful discussions. M.W., L.M. and M.B.
are supported by the Villum Foundation (Research Grant No. 25310). This
project has received funding from the European Union’s Horizon 2020 research
and innovation program under the Marie Skłodowska-Curie Grant Agreement No.
847523 “INTERACTIONS.” C.-M.C. acknowledges the support by the Ministry of
Science and Technology (MOST) under Grant No. 111-2112-M-110-006-MY3, and by
the Yushan Young Scholar Program under the Ministry of Education (MOE) in
Taiwan. This work was supported by Q@TN, the joint lab between University of
Trento, FBK—Fondazione Bruno Kessler, INFN—National Institute for Nuclear
Physics, and CNR—National Research Council.

##  Appendix A Spin relaxation: comparison with NRG

Figure 5: Comparison between the spin relaxation dynamics obtained with the
method presented in this paper (solid lines) and the data extracted from
Fig.4b of Ref. [25] (empty circles). data corresponding to the same ratio
$U/\Gamma$ have the same color.

As a further benchmark for our code, in this appendix we compare explicitly
the results obtained with our MPS method with the seminal time-dependent NRG
(TD-NRG) approach introduced in by Anders and Schiller in Ref. [25]. Here we
fix the impurity-lead coupling $\Gamma=0.08t_{0}$ and change the interaction
strength $U$, while the chemical potential is set at the particle-hole
symmetric point. The leads have a length $L=100$ and a decay length $\xi=16$.
As usual, we use a SVD truncation error of $5\cdot 10^{-8}$ and a maximum bond
dimension of 2000.

In Fig. 5 we show the spin relaxation of the impurity, initially prepared in
the $|\uparrow\rangle$ state, in a zero-bias quench. Curves corresponding to
different values of $U/\Gamma$ nicely collapse on top of each other if the
time is measured in units of the Kondo time scale ${\sf t}_{K}$ and we get
good agreement with the data extracted from Fig.4b of Ref. [25] using the
software WebPlotDigitizer [61]. Notice that we have rescaled their y-axis by a
factor of 4 to match our definition of the magnetization. The main differences
between the two approaches occur at small and large times; since we use the
same time step $\Delta t=0.5\hbar/t_{0}$ in the TDVP algorithm for all values
of $U/\Gamma$, we don’t have a good resolution at ${\sf t}\ll{\sf t}_{K}$, in
particular when the Kondo time is small. Hence, our simulations do not capture
precisely the initial plateau of the magnetization dynamics, as it requires a
smaller time step. At longer times, instead, we observe a relaxation of the
magnetization to smaller values, even though the entanglement growth prevents
us to reach timescale greater than $\sim 10{\sf t}_{K}$. This discrepancy
might be due to the different ways finite-size effects affect TD-NRG and TDVP
algorithms.

## References

  * Goldhaber-Gordon  _et al._ [1998] D. Goldhaber-Gordon, H. Shtrikman, D. Mahalu, D. Abusch-Magder, U. Meirav, and M. A. Kastner, Kondo effect in a single-electron transistor, Nature 391, 156 (1998).
  * Cronenwett  _et al._ [1998] S. M. Cronenwett, T. H. Oosterkamp, and L. P. Kouwenhoven, A tunable kondo effect in quantum dots, Science 281, 540 (1998).
  * Nygård  _et al._ [2000] J. Nygård, D. H. Cobden, and P. E. Lindelof, Kondo physics in carbon nanotubes, Nature 408, 342 (2000).
  * Jespersen  _et al._ [2006] T. S. Jespersen, M. Aagesen, C. Sørensen, P. E. Lindelof, and J. Nygård, Kondo physics in tunable semiconductor nanowire quantum dots, Phys. Rev. B 74, 233304 (2006).
  * Riegger  _et al._ [2018] L. Riegger, N. Darkwah Oppong, M. Höfer, D. R. Fernandes, I. Bloch, and S. Fölling, Localized magnetic moments with tunable spin exchange in a gas of ultracold fermions, Phys. Rev. Lett. 120, 143601 (2018).
  * Kanász-Nagy  _et al._ [2018] M. Kanász-Nagy, Y. Ashida, T. Shi, C. Moca, T. N. Ikeda, S. Fölling, J. I. Cirac, G. Zaránd, and E. A. Demler, Exploring the anisotropic kondo model in and out of equilibrium with alkaline-earth atoms, Phys. Rev. B 97, 155156 (2018).
  * Haldane [1978] F. D. M. Haldane, Theory of the atomic limit of the Anderson model. I. Perturbation expansions re-examined, J. Phys. C: Solid State Phys. 11, 5015 (1978).
  * Werner  _et al._ [2009] P. Werner, T. Oka, and A. J. Millis, Diagrammatic monte carlo simulation of nonequilibrium systems, Phys. Rev. B 79, 035320 (2009).
  * Werner  _et al._ [2010] P. Werner, T. Oka, M. Eckstein, and A. J. Millis, Weak-coupling quantum monte carlo calculations on the keldysh contour: Theory and application to the current-voltage characteristics of the anderson model, Phys. Rev. B 81, 035108 (2010).
  * Gull  _et al._ [2011] E. Gull, D. R. Reichman, and A. J. Millis, Numerically exact long-time behavior of nonequilibrium quantum impurity models, Phys. Rev. B 84, 085134 (2011).
  * Profumo  _et al._ [2015] R. E. V. Profumo, C. Groth, L. Messio, O. Parcollet, and X. Waintal, Quantum monte carlo for correlated out-of-equilibrium nanoelectronic devices, Phys. Rev. B 91, 245154 (2015).
  * Krivenko  _et al._ [2019] I. Krivenko, J. Kleinhenz, G. Cohen, and E. Gull, Dynamics of kondo voltage splitting after a quantum quench, Phys. Rev. B 100, 201104 (2019).
  * Kleinhenz  _et al._ [2022] J. Kleinhenz, I. Krivenko, G. Cohen, and E. Gull, Kondo cloud in a one-dimensional nanowire, Phys. Rev. B 105, 085126 (2022).
  * Anders [2008] F. B. Anders, Steady-state currents through nanodevices: A scattering-states numerical renormalization-group approach to open quantum systems, Phys. Rev. Lett. 101, 066804 (2008).
  * Souto  _et al._ [2018] R. S. Souto, R. Avriller, A. L. Yeyati, and A. Martín-Rodero, Transient dynamics in interacting nanojunctions within self-consistent perturbation theory, New Journal of Physics 20, 083039 (2018).
  * Al-Hassanieh  _et al._ [2006] K. A. Al-Hassanieh, A. E. Feiguin, J. A. Riera, C. A. Büsser, and E. Dagotto, Adaptive time-dependent density-matrix renormalization-group technique for calculating the conductance of strongly correlated nanostructures, Phys. Rev. B 73, 195304 (2006).
  * Dias da Silva  _et al._ [2008] L. G. G. V. Dias da Silva, F. Heidrich-Meisner, A. E. Feiguin, C. A. Büsser, G. B. Martins, E. V. Anda, and E. Dagotto, Transport properties and Kondo correlations in nanostructures: Time-dependent DMRG method applied to quantum dots coupled to Wilson chains, Phys. Rev. B 78, 195317 (2008).
  * Heidrich-Meisner  _et al._ [2009] F. Heidrich-Meisner, A. E. Feiguin, and E. Dagotto, Real-time simulations of nonequilibrium transport in the single-impurity Anderson model, Phys. Rev. B 79, 235336 (2009).
  * Schwarz  _et al._ [2018] F. Schwarz, I. Weymann, J. von Delft, and A. Weichselbaum, Nonequilibrium steady-state transport in quantum impurity models: A thermofield and quantum quench approach using matrix product states, Phys. Rev. Lett. 121, 137702 (2018).
  * Manaparambil  _et al._ [2022] A. Manaparambil, A. Weichselbaum, J. von Delft, and I. Weymann, Nonequilibrium spintronic transport through kondo impurities, Phys. Rev. B 106, 125413 (2022).
  * Kohn and Santoro [2022] L. Kohn and G. E. Santoro, Quench dynamics of the Anderson impurity model at finite temperature using matrix product states: entanglement and bath dynamics, Journal of Statistical Mechanics: Theory and Experiment 2022, 063102 (2022).
  * Thoenniss  _et al._ [2023a] J. Thoenniss, M. Sonner, A. Lerose, and D. A. Abanin, Efficient method for quantum impurity problems out of equilibrium, Phys. Rev. B 107, L201115 (2023a).
  * Eckel  _et al._ [2010] J. Eckel, F. Heidrich-Meisner, S. G. Jakobs, M. Thorwart, M. Pletyukhov, and R. Egger, Comparative study of theoretical methods for non-equilibrium quantum transport, New J. Phys. 12, 043042 (2010).
  * Chung  _et al._ [2022] C.-M. Chung, M. M. Wauters, and M. Burrello, Matrix product state simulations of quantum quenches and transport in coulomb blockaded superconducting devices, Phys. Rev. B 106, 094308 (2022).
  * Anders and Schiller [2005] F. B. Anders and A. Schiller, Real-time dynamics in quantum-impurity systems: A time-dependent numerical renormalization-group approach, Phys. Rev. Lett. 95, 196801 (2005).
  * Aoki  _et al._ [2014] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Nonequilibrium dynamical mean-field theory and its applications, Rev. Mod. Phys. 86, 779 (2014).
  * Weichselbaum  _et al._ [2009] A. Weichselbaum, F. Verstraete, U. Schollwöck, J. I. Cirac, and J. von Delft, Variational matrix-product-state approach to quantum impurity models, Phys. Rev. B 80, 165117 (2009).
  * Chien  _et al._ [2014] C.-C. Chien, M. Di Ventra, and M. Zwolak, Landauer, Kubo, and microcanonical approaches to quantum transport and noise: A comparison and implications for cold-atom dynamics, Phys. Rev. A 90, 023624 (2014).
  * Bertini  _et al._ [2016] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Transport in out-of-equilibrium $xxz$ chains: Exact profiles of charges and currents, Phys. Rev. Lett. 117, 207201 (2016).
  * Essler and Fagotti [2016] F. H. L. Essler and M. Fagotti, Quench dynamics and relaxation in isolated integrable quantum spin chains, Journal of Statistical Mechanics: Theory and Experiment 2016, 064002 (2016).
  * Lieb and Robinson [1972] E. H. Lieb and D. W. Robinson, The finite group velocity of quantum spin systems, Commun. Math. Phys. 28, 251 (1972).
  * Bonnes  _et al._ [2014] L. Bonnes, F. H. L. Essler, and A. M. Läuchli, “light-cone” dynamics after quantum quenches in spin chains, Phys. Rev. Lett. 113, 187203 (2014).
  * Viti  _et al._ [2016] J. Viti, J.-M. Stéphan, J. Dubail, and M. Haque, Inhomogeneous quenches in a free fermionic chain: Exact results, EPL (Europhysics Letters) 115, 40011 (2016).
  * Ljubotina  _et al._ [2019] M. Ljubotina, S. Sotiriadis, and T. Prosen, Non-equilibrium quantum transport in presence of a defect: the non-interacting case, SciPost Phys. 6, 4 (2019).
  * Boulat  _et al._ [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).
  * Güttge  _et al._ [2013] F. Güttge, F. B. Anders, U. Schollwöck, E. Eidelstein, and A. Schiller, Hybrid NRG-DMRG approach to real-time dynamics of quantum impurity systems, Phys. Rev. B 87, 115115 (2013).
  * He and Millis [2017] Z. He and A. J. Millis, Entanglement entropy and computational complexity of the Anderson impurity model out of equilibrium: Quench dynamics, Phys. Rev. B 96, 085107 (2017).
  * Thoenniss  _et al._ [2023b] J. Thoenniss, A. Lerose, and D. A. Abanin, Nonequilibrium quantum impurity problems via matrix-product states in the temporal domain, Phys. Rev. B 107, 195101 (2023b).
  * Ng  _et al._ [2023] N. Ng, G. Park, A. J. Millis, G. K.-L. Chan, and D. R. Reichman, Real-time evolution of anderson impurity models via tensor network influence functionals, Phys. Rev. B 107, 125103 (2023).
  * Rosch [2012] A. Rosch, Wilson chains are not thermal reservoirs, Eur. Phys. J. B 85, 1 (2012).
  * Haegeman  _et al._ [2011] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Time-dependent variational principle for quantum lattices, Phys. Rev. Lett. 107, 070601 (2011).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Unifying time evolution and optimization with matrix product states, Phys. Rev. B 94, 165116 (2016).
  * Rams and Zwolak [2020] M. M. Rams and M. Zwolak, Breaking the entanglement barrier: Tensor network simulation of quantum transport, Phys. Rev. Lett. 124, 137701 (2020).
  * Wauters  _et al._ [2023b] M. M. Wauters, C.-M. Chung, L. Maffi, and M. Burrello, Topological kondo model out of equilibrium, Phys. Rev. B 108, L220302 (2023b).
  * Fishman  _et al._ [2022] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022).
  * Note [1] The source code can be found in the repository https://github.com/chiamin/QuenchTranportTwoChains.
  * Lobaskin and Kehrein [2005] D. Lobaskin and S. Kehrein, Crossover from nonequilibrium to equilibrium behavior in the time-dependent kondo model, Phys. Rev. B 71, 193303 (2005).
  * Cavalcante  _et al._ [2023] M. F. Cavalcante, R. G. Pereira, and M. C. O. Aguiar, Quench dynamics of the kondo effect: Transport across an impurity coupled to interacting wires, Phys. Rev. B 107, 075110 (2023).
  * Wang  _et al._ [2008] X. Wang, C. D. Spataru, M. S. Hybertsen, and A. J. Millis, Electronic correlation in nanoscale junctions: Comparison of the gw approximation to a numerically exact solution of the single-impurity anderson model, Phys. Rev. B 77, 045119 (2008).
  * Kaminski  _et al._ [2000] A. Kaminski, Yu. V. Nazarov, and L. I. Glazman, Universality of the Kondo effect in a quantum dot out of equilibrium, Phys. Rev. B 62, 8154 (2000).
  * Zitko and Pruschke [2009] R. Zitko and T. Pruschke, Energy resolution and discretization artifacts in the numerical renormalization group, Phys. Rev. B 79, 085106 (2009).
  * Zitko [2021] R. Zitko, Nrg ljubljana (2021).
  * Béri and Cooper [2012] B. Béri and N. R. Cooper, Topological Kondo effect with Majorana fermions, Phys. Rev. Lett. 109, 156803 (2012).
  * Béri [2013] B. Béri, Majorana-Klein hybridization in topological superconductor junctions, Phys. Rev. Lett. 110, 216803 (2013).
  * Altland and Egger [2013] A. Altland and R. Egger, Multiterminal Coulomb-Majorana junction, Phys. Rev. Lett. 110, 196401 (2013).
  * Kanne  _et al._ [2022] T. Kanne, D. Olsteins, M. Marnauza, A. Vekris, J. C. Estrada Saldaña, S. Loric, R. D. Schlosser, D. Ross, S. Csonka, K. Grove-Rasmussen, and J. Nygård, Double nanowires for hybrid quantum devices, Advanced Functional Materials 32, 2107926 (2022).
  * Vekris  _et al._ [2022] A. Vekris, J. C. Estrada Saldaña, T. Kanne, T. Hvid-Olsen, M. Marnauza, D. Olsteins, M. M. Wauters, M. Burrello, J. Nygård, and K. Grove-Rasmussen, Electronic transport in double-nanowire superconducting islands with multiple terminals, Nano Letters 22, 5765 (2022), pMID: 35833741, https://doi.org/10.1021/acs.nanolett.2c01161 .
  * Souto  _et al._ [2022] R. S. Souto, M. M. Wauters, K. Flensberg, M. Leijnse, and M. Burrello, Multiterminal transport spectroscopy of subgap states in coulomb-blockaded superconductors, Phys. Rev. B 106, 235425 (2022).
  * Oreg and von Oppen [2020] Y. Oreg and F. von Oppen, Majorana zero modes in networks of cooper-pair boxes: Topologically ordered states and topological quantum computation, Annual Review of Condensed Matter Physics 11, 397 (2020).
  * Beenakker [2020] C. W. J. Beenakker, Search for non-Abelian Majorana braiding statistics in superconductors, SciPost Phys. Lect. Notes , 15 (2020).
  * Rohatgi [2022] A. Rohatgi, Webplotdigitizer: Version 4.6 (2022).