# Out-of-equilibrium dynamics in a quantum impurity model: Numerics for particle transport and entanglement entropy

###### Abstract

We investigate the out-of-equilibrium properties of a simple quantum impurity model, the interacting resonant level model. We focus on the scaling regime, where the bandwidth of the fermions in the leads is larger than all the other energies, so that the lattice and the continuum versions of the model become equivalent. Using time-dependent density matrix renormalization group simulations initialized with states having different densities in the two leads we extend the results of Boulat, Saleur and Schmitteckert [Phys. Rev. Lett. 101, 140601 (2008)] concerning the current-voltage (I-V) curves, for several values of the interaction strength U. We estimate numerically the Kondo scale T_{B} and the exponent b(U) associated to the tunneling of the fermions from the leads to the dot. Next we analyze the quantum entanglement properties of the steady states. We focus in particular on the entropy rate \alpha, describing the linear growth with time of the bipartite entanglement in the system. We show that, as for the current, \alpha/T_{B} is described by some function of U and of the rescaled bias V/T_{B}. Finally, the spatial structure of the entropy profiles is discussed.

## I Introduction

The dynamical properties of out-of-equilibrium quantum many-body systems are a major topic in condensed-matter physics Polkovnikov et al. (2011); Eisert, Friesdorf, and Gogolin (2011). Although the equilibrium properties of one-dimensional (1d) interacting problems are well understood in many cases, thanks, for instance, to some theoretical techniques such as the Bethe Ansatz, conformal field theory or matrix-product states (MPS) numerics Schollwöck (2011), the physics that takes place in out-of-equilibrium situations represents very rich domain, and is much less understood.

Quantum impurity problems are among the simplest quantum many-body systems, but they nevertheless harbor many interesting phenomena, many open questions, and they represent a very useful playground to develop new theoretical ideas and methods. They are also – of course – relevant to describe many experimental situations, from the Kondo effect in metals Hewson (1993) to transport in nanostructures such as quantum dots or point contacts.

In this work we consider a well known impurity model - the interacting resonant level model (IRLM) Wiegmann and Finkel’shtein (1978). The model describes two semi-infinite leads with spinless fermions that are coupled to some resonant level (dot) via some tunneling and some interaction [see Eq. (4) below]. We are interested here in the transport properties of the system, and wish to describe quantitatively the steady states that appear when some particle current is flowing through the dot. How do we address these questions, without relying on the linear response theory, nor using some perturbative scheme that would assume that the system is close to equilibrium and/or that the interactions are weak ? Thanks to the integrability of the IRLM Filyov and Wiegmann (1980), several remarkable exact results have been obtained concerning non-equilibrium steady states, such as the current-voltage (I-V) characteristic for some special (so-called “self-dual”) value of the interaction strength Boulat, Saleur, and Schmitteckert (2008) (see also Fendley, Ludwig, and Saleur (1995a, b); Mehta and Andrei (2006)). However, apart from this special point and the noninteracting case, simple quantities such as the I-V curve are not known analytically. From this point of view numerical simulations are invaluable and complementary to the analytical approaches.

This study indeed aims at providing accurate numerical data concerning the so-called scaling regime of the lattice model (i.e. all the energies are small compared to the bandwidth in the leads), where many quantities become universal and can be quantitatively compared to the field theory results (continuum limit).

To simulate transport, a useful setup is to prepare a large but finite isolated system in an initial state, at t=0, where two spatial regions – say the left and the right leads – have different particle densities. Starting from such an inhomogeneous state, the Hamiltonian evolution of the wave function will put the particles in motion. This first leads to some transient regime where some current starts flowing from one side to the other. For times that are long compared to the microscopic time scales, but smaller than the time required for an elementary excitation to propagate through the whole system, we expect on physical grounds that some quasi-steady states should be realized. This type of approach, where one follows numerically the real time evolution starting from an initial state with finite density bias, has already been used to investigate several interacting 1d systems, like XXZ spin-\frac{1}{2} chains Gobert et al. (2005); Sabetta and Misguich (2013), or impurity models Schmitteckert (2004); Schneider and Schmitteckert (2006); Boulat, Saleur, and Schmitteckert (2008). These simulations have been performed using the time-evolving block decimation and time-dependent density matrix renormalization group (DMRG) Vidal (2004); White and Feiguin (2004); Daley et al. (2004), where the many-body wave function of the system is encoded as a matrix-product state.

Our objective is first to refine and extend some of the previous numerical studies concerning the particle current in the IRLM Boulat, Saleur, and Schmitteckert (2008), and then to focus on the entanglement entropy. Even though there have been several studies on the entanglement properties of quantum impurity models, most of these works focused on the ground state (for a review see Laflorencie (2016)). Here we will analyze in detail the linear growth of the entanglement entropy with time, and its spatial structure in the steady regime.

The plan of the paper is as follows. Section. II presents the model, the initial state, and describes qualitatively the evolution of three quantities of interest: particle density, particle current, and entanglement (von Neumann) entropy. The central results are then presented in Sec. III. This section is devoted to the steady state, and presents some quantitative analysis of the numerics for the steady current and entropy rate, focusing on the scaling regime of the model. In particular we confirm (Sec. III.2), following Ref. Boulat, Saleur, and Schmitteckert (2008), that in this limit the current I is some universal function of V/T_{B} and U, where V is the initial bias, U the interaction strength, and T_{B} is Kondo crossover scale that we evaluate numerically. In Sec. III.3 we present a similar analysis for the scaling of the entropy rate, and we show that it can also be analyzed in term of some universal functions of V/T_{B} that we compute numerically for several values of U. The shape of the out-of-equilibrium entanglement profile is discussed in Sec. III.4.

## II Model and time evolution

### II.1 Hamiltonian

We consider a lattice version of the IRLM (Fig. 1), which can be defined as

\displaystyle H_{\rm IRLM} | \displaystyle=H_{A}+H_{B}+H_{d}, | (1) | |||

\displaystyle H_{A} | \displaystyle=-J\sum_{r=-N/2}^{-2}\left(c^{\dagger}_{r}c_{r+1}+\text{H.c}% \right), | (2) | |||

\displaystyle H_{B} | \displaystyle=-J\sum_{r=1}^{N/2-1}\left(c^{\dagger}_{r}c_{r+1}+\text{H.c}% \right), | (4) | |||

\displaystyle H_{d} | \displaystyle=-J^{\prime}\sum_{r=-1}^{0}\left(c^{\dagger}_{r}c_{r+1}+\text{H.c% }\right)+ | ||||

\displaystyle+U\sum_{r=\pm 1}\left(c^{\dagger}_{r}c_{r}-\frac{1}{2}\right)% \left(c^{\dagger}_{0}c_{0}-\frac{1}{2}\right), |

where H_{A/B} describes the kinetic energy of free spinless fermions in the left and right leads, and H_{d} models the tunneling from the leads to the dot (site at r=0), as well as the density-density interaction between the dot and the ends of the leads (r=-1 and 1).

As discussed later we will focus on the regime where the hopping amplitude J^{\prime} (or tunneling strength) between the leads and the dot (r=0) is much smaller than the bandwidth 4J of the kinetic energy in the leads, i.e J^{\prime}\ll J. From now we take J=1=\hbar, thus defining the unit of time and energy.

### II.2 Initial state

We choose the initial state |\psi(t=0)\rangle to be the ground state of H_{0}=H_{\rm IRLM}+H_{\rm bias}, where H_{\rm bias} is an inhomogeneous chemical potential (or voltage) that induces different densities on the left and on the right leads:

H_{\rm bias}=\frac{1}{2}V\sum_{r=-N/2}^{N/2-1}\tanh(r/w)c^{\dagger}_{r}c_{r}. | (5) |

In the left (right) lead the chemical potential is thus equal to V/2 (-V/2) sufficiently far from the dot. This bias induces different initial densities \langle c^{\dagger}_{r}c_{r}\rangle=\frac{1}{2}\pm m_{0} in the bulk of the leads at t=0 (blue horizontal line in Fig. 2). For an infinite system, N\to\infty, the Fermi momenta k_{F}^{+} (k_{F}^{-}) in the left (right) lead is set by 2J\cos(k_{F}^{\pm})=\pm V/2, and these are related to the density difference m_{0} through k_{F}^{\pm}=\pi(\pm m_{0}+1/2).

As done in previous studies Boulat, Saleur, and Schmitteckert (2008), the voltage drop in Eq. (5) is spatially smeared over \sim w sites in the vicinity of the dot. This has the effect of producing an initial state with smoother density in the vicinity of r=0 and turns out to accelerate the convergence to a steady state. We typically use w=10. For the same reason, the initial state is prepared with J^{\prime}=J, that is uniform hopping amplitudes throughout the chain. In addition, H_{0} is chosen to be noninteracting (U=0), and the interactions are switched on for t>0.

### II.3 Unitary evolution at t>0

For t>0 the wave function evolves using H_{\rm IRLM}, with 0<J^{\prime}<1, and without the voltage bias term. The calculations are performed using a time-dependent DMRG algorithm, implemented using the C++ iTensor library ite (n 20). The evolution operator U=\exp(-i\tau H) for a time step \tau is approximated by a matrix-product operator (MPO) Zaletel et al. (2015), using a fourth-order Trotter scheme [Eqs.(10) and (LABEL:eq:trotter4)] and \tau=0.2 (unless specified otherwise). The system sizes are of the order of 200 sites and the largest times are of the order of t\simeq 100. Additional details about the simulations are given in Appendix A.

### II.4 Particle density

The particle density is defined by \rho(r,t)=\langle\psi(t)|c^{\dagger}_{r}c_{r}|\psi(t)\rangle, and we can equivalently use a spin language, with the magnetization S^{z}(r,t)=\rho(r,t)-\frac{1}{2}. A typical evolution of the density profile is shown in Fig. 2. It shows how the initial profile at t=0 gives rise to two propagating fronts (one to the left and one to the right), forming a “light cone”, and how some steady region form in the center. When the time is large enough, two regions with quasi homogeneous densities develop on both sides of the dot. The densities in the steady regions of the lead can be written as \rho=\frac{1}{2}\pm m, and the density difference m between both sides of the dot can be computed exactly in the free-fermion case U=0. The result reads as:

m=\int_{k_{F}^{-}}^{k_{F}^{+}}\frac{dk}{2\pi}\mathcal{R}(\epsilon(k)) | (6) |

where \mathcal{R}(\epsilon) is the reflexion coefficient for an incident particle with energy \epsilon (more details in Appendix. B.2).

### II.5 Particle current

On a given bond the expectation value of the current operator is I(r,t)=2J_{r}{\rm Im}\langle\psi(t)|c^{\dagger}_{r}c_{r+1}|\psi(t)\rangle, where J_{r} is the hopping amplitude between the sites located at r and r+1. We thus have J_{r}=J^{\prime} for r=-1 and r=0 (hopping to the dot), and J_{r}=J=1 otherwise (in the leads). This definition ensures the proper charge conservation equation, \frac{d}{dt}\langle\psi(t)|c^{\dagger}_{r}c_{r}|\psi(t)\rangle=I(r-1,t)-I(r,t). When no position r is given, I(t) refers to the averaged current on both sites of the dot: I(t)=\frac{1}{2}\left(I(-1,t)+I(0,t)\right).

A typical evolution of the current profile is shown in the upper panel of Fig. 3. As for the density, two propagating fronts are visible. In the center of the system one observes, at sufficiently long times, the emergence of a spatial region where the current is almost constant in space and time. This is where some steady value of the current can be defined. Interestingly, the current in the (left-moving and right-moving) front regions reaches values that are significantly larger than in the steady regions.

The way the current in the center of the system approaches its steady value, after some damped oscillations, is shown in the upper panel of Fig. 4.

In many cases the oscillations which appear in the transient regime have a well-defined period T_{\rm osc}. It was argued in Ref. Schneider and Schmitteckert (2006) that this period is simply given by the bias: T_{\rm osc}=4\pi/V. Our data are in agreement with this result, from U=0 to large values of U.

### II.6 Entanglement entropy

We denote by S(t,R+\frac{1}{2}) the von Neumann entanglement entropy of the set A of sites located to the left of site R, that is A=\{-N/2,\dots,R\}. When no position is specified, S(t) refers to the entanglement entropy S(t,-\frac{1}{2}) of the entire left (or right) lead.

The lower panel of Fig. 3 illustrates how the entanglement profile evolves. The most striking feature is the rapid growth of the entanglement entropy, which turns out to be linear in time for a given r inside the “light cone” (see the lower panel in Fig. 4). This linear growth of the entropy is well known in the situations where some steady current is flowing through an impurity (or defect). One can in particular mention the analogous case of a weak bond connecting two free leads, studied in detail in Ref. Eisler and Peschel (2012). In such situations, a quantity of interest is the entropy rate, defined as \alpha=\frac{d}{dt}S(t).

Since the computational cost of MPS-based methods grows exponentially with the amount of bipartite entanglement entropy in the system, this linear entropy growth severely limits the longest times we can reach in the simulations. This should be contrasted, for instance, with the slower logarithmic growth of the entropy in the case where the two leads are connected to the dot in the absence of any bias Kennes, Meden, and Vasseur (2014); Vasseur and Saleur (2017). A logarithmic growth is in fact generic for local quenches in critical one dimensional systems Calabrese and Cardy (2007); Stéphan and Dubail (2011); Eisler and Peschel (2012). The entropy growth is also logarithmic in the case of an XXZ spin chain (|\Delta|<1), with a translation invariant Hamiltonian, which is set in a current-carrying state using some domain-wall initial condition (equivalent to some non zero bias) Gobert et al. (2005); Sabetta and Misguich (2013).

The physical origin of the finite rate \alpha is easy to understand in a noninteracting and semiclassical picture. Each incoming particle at energy \epsilon has a finite probability \mathcal{T}(\epsilon) to be transmitted to the other side of the dot, and a probability \mathcal{R}(\epsilon)=1-\mathcal{T}(\epsilon) to be reflected (more details in Appendix B.3). After such a scattering event, the wave packet of the particle is split in two parts, one on each side of the dot, propagating in opposite directions. In other words, the state of this particle is a quantum superposition of two terms, one in which the particle is in the left lead, and another one where the particle is in the right lead. So, each such event contributes by an amount \delta S=-\mathcal{T}(\epsilon)\ln\mathcal{T}(\epsilon)-\mathcal{R}(\epsilon)% \ln\mathcal{R}(\epsilon) to the entanglement entropy between the two leads. In presence of a finite steady current there is finite charge transmitted per unit of time, and hence a linear growth of the entropy S(t)\sim\alpha t [except if \mathcal{R}(\epsilon)=0, as for J^{\prime}=J]. In such a picture, the entanglement is directly related to quantum fluctuations of the transmitted charge, present as soon as \mathcal{T}(\epsilon) is different from 1 and from 0, and this is nothing but a manifestation of the relation between entanglement and charge fluctuations in free particle systems Klich and Levitov (2009); Song et al. (2012, 2011). In contrast, if U\neq 0, the entanglement growth is a priori not directly related to the partial transfer of the particles. We will indeed see in Sec. III.3 that one can be in a situation where I goes to zero while the rate \alpha stays finite.

The semiclassical description at U=0 can also be used to understand qualitatively the triangular shape of the entropy profile. Indeed,
in the limit where the bias V is small compared to the bandwidth, all the relevant excitations propagate at the same group velocity (\pm v_{F}=\pm 2J). In that case, the degrees of freedom which contribute to the entanglement entropy form a “train” of left-moving wave packets with momenta centered around -k_{F}, and another train with right-moving wave packets centered at +k_{F} ^{1}^{1}1The average spacing \lambda between the packets is proportional to the inverse of their momentum width: \lambda=2\pi(\Delta k)^{-1} . The momentum width is related
to the energy width, v_{F}\Delta k=\Delta\epsilon, and this energy width is nothing but the bias, \Delta\epsilon=V.
We thus have v_{F}\Delta k=V. The number of incident particles per unit of time
is thus v_{F}/\lambda=V/(2\pi), and this is consistent with the small V limit of the current given in Eq. 20.
Note that the actual distribution is in fact Poissonian.
.
The important point is that each left-moving particle is entangled with one right-moving partner, located on the other side of the dot, at the same distance.
When performing a partition of the chain at a given time t and at a given position R, the entanglement entropy that is predicted by the semiclassical description simply depends on the number
of entangled pairs which are separated by the partition. It is then straightforward to see that this leads to a triangular entropy profile,
with a spatial extension ranging from R=-v_{F}t to +v_{F}t, and a height equal to \alpha t, with the rate \alpha given in Eq. (26). Although this classical picture with noninteracting particles
does not apply to the interacting case, the numerical simulations show that the triangular shape of the entropy profiles is a robust feature, at least far enough from the dot.
The integrability of the IRLM implies that, in some sense, a particle picture is still applicable, even in presence of strong interactions. This property was crucial
to derive the results in Refs. Boulat, Saleur, and Schmitteckert (2008); Fendley, Ludwig, and Saleur (1995a, b), and it might be related to the triangular profile observed here.
Closer to the dot, the profile is, however, not triangular, and this will be analyzed in Sec. III.4.

It should finally be noted that this linear entropy growth is what makes this type of simulation difficult, since it forces the matrix dimensions in the MPS representation of the wave function to grow exponentially with time. Some more details on this point can be found in Appendix A.2.

## III Nonequilibrium steady states

We discuss here the properties of the steady region which grows in the center of the system, and where local observables asymptotically become independent of time.

### III.1 Reminder on the scaling regime and the continuum limit of the IRLM

As mentioned in the Introduction, we focus on the regime where
the free-fermion bandwidth W=4J is larger than all the other energies in the problem, namely, J^{\prime}, V, and U ^{2}^{2}2In practice we take J=1 and restrict our numerical simulations to V\lesssim 2, J^{\prime}\lesssim 0.3 and U\leq 2.5. Finite-J^{\prime} corrections start to be visible for J^{\prime}=0.5. As for U, our data at U=4 are clearly outside the scaling regime..

In this regime, the microscopic details of the leads (like the band structure) are irrelevant, except for the Fermi velocity (v_{F}=2J), and their gapless low-energy excitations are described in the continuum limit in terms of simple scale-invariant Hamiltonians for left- and right-moving relativistic fermions:

\displaystyle H_{A}^{c} | \displaystyle= | \displaystyle iv_{F}\int_{-\infty}^{0}dr\left[\psi_{L}^{\dagger}(r)\partial_{r% }\psi_{L}(r)-\psi_{R}^{\dagger}(r)\partial_{r}\psi_{R}(r)\right] | (7) | ||

\displaystyle H_{B}^{c} | \displaystyle= | \displaystyle iv_{F}\int_{0}^{\infty}dr\left[\psi_{L}^{\dagger}(r)\partial_{r}% \psi_{L}(r)-\psi_{R}^{\dagger}(r)\partial_{r}\psi_{R}(r)\right]. | (8) |

H_{A}^{c} describes the continuum limit of the left lead (r<0), H_{B}^{c} describes that of the right lead (r>0), and \psi_{L} and \psi_{R} are the left- and right-moving fermion annihilation operators. The coupling to the dot then takes the form

\displaystyle H_{d}^{c} | \displaystyle= | \displaystyle-J_{c}^{\prime}\left(\psi_{L}(0)+\psi_{R}(0)\right)d^{\dagger}+{% \rm H.c.} | (9) | ||

\displaystyle+U_{c}\left(:\psi_{L}(0)^{\dagger}\psi_{L}(0):+:\psi_{R}(0)^{% \dagger}\psi_{R}(0):\right) | |||||

\displaystyle\times\left(d^{\dagger}d-\frac{1}{2}\right), |

where d is the fermion operator associated to the dot. The analysis of this model then usually proceeds by “unfolding” the two semi-infinite leads, giving two infinite right-moving Fermi wires, but we will not pursue this here.

Since the tunneling to the dot is a relevant interaction, a Kondo energy scale T_{B} appears when the leads are connected to the dot, and most quantities are expected to follow some single-parameter scaling with T_{B}. At energies that are small compared to the crossover scale T_{B}, the dot is hybridized with the leads, and the system effectively appears as a single chain (so-called “healing”), whereas the wires appear to be almost disconnected at energies much larger than T_{B}. As for the interaction U, it corresponds to a marginal operator and therefore changes continuously the critical properties of the model. Although this is well established for equilibrium properties, it is less obvious that T_{B} also rules the out-of-equilibrium properties. Such behavior is nevertheless verified for the current I, which, for a given U, takes the form I/T_{B}=f(V/T_{B}) Boulat, Saleur, and Schmitteckert (2008). The role played by T_{B} in the dynamics of the IRLM has also been investigated in the absence of any bias (V=0), in a local quench setup when the leads are abruptly connected to the dot at t=0 Kennes, Meden, and Vasseur (2014); Vasseur and Saleur (2017).

In Sec. II.6 we show that the stationary rate \alpha=\frac{d}{dt}S at which the entanglement entropy in the center grows with time obeys a similar scaling form.

The energy T_{B} is known to scale as \sim J^{\prime\frac{1}{1-h}}, where h is the scaling dimension of the operator which, in the continuum description, describes the tunneling from the leads to the dot. This dimension h depends on the interaction through h=\frac{1}{4}+\left(\frac{1}{2}-\frac{U_{c}}{2\pi}\right)^{2}, where U_{c} is the interaction strength in the continuum limit Boulat, Saleur, and Schmitteckert (2008). Finding the precise way U_{c} depends on the microscopic parameter U of the lattice model would require to follow exactly the renormalization-group flow going from the lattice model to the infra red fixed point, but there is no known method to do this exactly. The analytical result of Ref. Boulat, Saleur, and Schmitteckert (2008) was obtained for the special value U_{c}=\pi, where the model has some self-duality property and can be mapped onto the boundary sine-Gordon model. Thanks to numerical simulations, it has been shown that the lattice model at U=2.0 has a continuum limit which is close to the self dual point, where the exponent h reaches a minimum Boulat, Saleur, and Schmitteckert (2008). Our data also confirm this result. We also improve quantitatively the connection between interaction strength U on the lattice, and its value U_{c} in the infrared limit.

The exponent h also appears in the limit of large (rescaled) bias V/T_{B}, where the steady current behaves as a power law: I\sim V^{-b} with b=1-2h Boulat, Saleur, and Schmitteckert (2008). This behavior will be checked numerically in detail in the next section (Sec. III.2). For a given U, this offers a simple way to extract h and b from the simulations, and then to define T_{B} for each value of J^{\prime}.

### III.2 Steady Current

Figure 5 shows the stationary current as a function of the bias, for a few values of U and J^{\prime}=0.08.
A log-log scale is used to visualize the power-law behavior of the current I\sim V^{-b(U)} at sufficiently large V/T_{B}, i.e small J^{\prime}.
The slope of the curve in log-log scale allows to determine the exponent b(U) and h(U)=\frac{1}{2}[1-b(U)]. The results of these fits are shown in Fig. 6
^{3}^{3}3The value J^{\prime}=0.08 appears to offer a good compromise to estimate b in our simulations.
Indeed, we need a small J^{\prime} to be in the scaling regime, but
the time to reach the steady regime (which increases when J^{\prime} decreases) should also not become too large compared to L/v_{F}..
Then the scale T_{B} is defined as T_{B}(J^{\prime},U)=(J^{\prime})^{\frac{1}{1-h(U)}}=(J^{\prime})^{\frac{2}{1+b%
(U)}} ^{4}^{4}4We adopt the convention where the numerical prefactor in the definition of T_{B} is set to 1.
. From the analysis of the model in the continuum, we know that the exponent b should reach a maximum value b=\frac{1}{2} (equivalent to h=\frac{1}{4}) at the self-dual point.
The maximum value we obtain (Fig. 6) is b=0.494, which gives an estimate on our precision on this quantity.

With the above T_{B} one can define a rescaled current I/T_{B} and rescaled bias V/T_{B}. As shown in Figs. 7 and 8, we then observe a relatively good collapse, on a single master curve, of the data sets corresponding to different J^{\prime} (for a given U). This indicates that the lattice model is indeed close to the scaling regime, characterized by a single energy scale T_{B}. For U=2.0 this has already been observed by Boulat et al. Boulat, Saleur, and Schmitteckert (2008), but thanks to longer simulations and larger systems, the present data have some higher precision and we could extend the I-V curves to larger values of V/T_{B} (beyond 100) and for several values of U from -0.1 to 3.

The fact that the current decreases with V at large bias for U>0, called negative differential conductance, is a remarkable phenomenon due to the interaction (for U\leq 0 the current is monotonically increasing), and has already been discussed in Boulat, Saleur, and Schmitteckert (2008).

For small values of U the exponent b describing the current suppression at large bias has been computed using some functional renormalization group method Karrasch et al. (2010) or with some self-consistent conserving approximation Vinkler-Aviv, Schiller, and Anders (2014). Using our convention for the lattice model, their result reads as b=2U/\pi+\mathcal{O}(U^{2}). As shown in Fig. 6, this is consistent with our simulations. Our results however appear to be slightly below this first order expansion in U. Together with the fact that the maximum of b is found to be slightly below 0.5, this indicates that our procedure slightly underestimates the exponent. This effect is presumably due to the fact that the calculations are performed using a finite J^{\prime} (0.08) and finite voltage V (up to \simeq 2).

### III.3 Entropy rate

We estimate the steady entropy rate \alpha, defined in Sec. II.6, by fitting the long-time part of the entanglement entropy data S(t)
^{5}^{5}5
Like for the current I(t), the entropy S(t) often shows some oscillatory part, \sim\cos(Vt/2+{\rm cst}) (see Fig. 15),
and we include such a term in the fitting function in order to improve the precision on \alpha..
As for the current, we consider the rescaled
entropy rate \alpha/T_{B} as a function of the rescaled bias V/T_{B}. The results, plotted in Figs. 9 and 10,
show that the data obtained for a given value U but for different values of J^{\prime} collapse quite well onto a single master curve.
The values of T_{B} used to construct the Figs. 9 and 10 are the same as those used to analyze the scaling of the current.
From this point of view, the quality of the collapse for the entropy rate
is quite remarkable since there is no adjustable parameter: the scale T_{B} was extracted from the large-bias behavior of current only, and for a single value of J^{\prime}=0.08. Note also that, to our knowledge, no exact result is known for the entropy rate when U\neq 0, even at the self-dual point.

The free-fermion result for the entropy rate (derived in Appendix B.3) converges to some finite constant, \alpha/T_{B}\to 2, at large rescaled bias. An important fact we learn from the Fig. 9 is that, in presence of interactions, \alpha/T_{B} also saturates to some finite value at large V/T_{B}. This limiting value appears to be smaller than 2 when U>0. From our data the large bias behavior of the entropy rate when U<0 is not simple to guess. It may diverge as V/T_{B}\to\infty, or it may saturate to some value larger than 2.

For U>0, while the current decreases to zero at large voltage (positive exponent b), the entropy rate stays finite. This observation is somehow counter intuitive since a vanishingly small amount of charge is transferred per unit of time from one wire to the other, while their entanglement entropy still grows at a finite rate. In this large bias regime it is possible that the entanglement is generated by the density-density interaction U, without any actual transfer of particles from one lead to the other. This question certainly deserves some further investigations.

### III.4 Stationary entropy profile

Several works have shown that, in Kondo-type problems, the entanglement entropy can be used to identify some spatial region of size \xi\sim T_{B}^{-1} around the impurity (sometimes called Kondo ”cloud”). See, for instance, Refs. Freton, Boulat, and Saleur (2013); Vasseur and Saleur (2017) concerning the IRLM, and Ref. Laflorencie (2016) for a more general review. While these studies investigated the entanglement in the ground state of the model, here we instead look at the entropy profile in some nonequilibrium steady state, in presence of a finite current.

Our results are summarized in Figs. 11 and Figs. 12. As discussed in Sec. II.6, the global shape of the profile is approximately triangular, with some maximum value that grows as \alpha t with a rate \alpha\sim\mathcal{O}(T_{B}). In order to analyze more precisely the long time limit of the profiles, we subtract the value of the maximum of each profile.

When plotted as a function of the rescaled distance r/\xi, with \xi=T_{B}^{-1}, the data for U=0 corresponding to different values of J^{\prime} (but same V/T_{B}) collapse onto a single curve, at least approximately. At distances from the dot which are large compared to \xi, the curve is linear. In this region the collapse is a consequence of the fact that the entropy rate \alpha, and thus the slope of the profile, scales as T_{B}. However, the curve shows a nontrivial structure for distances from the dot that are of the order of \xi, with some rounded maximum. We interpret this as some signature of the more complex correlations taking place in a nonequilibrium Kondo ’cloud’ of size T_{B}^{-1}.

The situation at U=2, displayed in Fig. 12 is more intriguing. The use of the rescaled distance r\cdot T_{B} still allows the data associated to different J^{\prime} (and thus different T_{B}) to collapse on a straight line for distances that are sufficiently large. But since the entanglement entropy is just a linear function of the distance to the dot when r is sufficiently large, this collapse only reflects the fact that the slope of the entropy profile is proportional to T_{B} (which we know already, since the entropy rate scales as T_{B} and the Fermi velocity is equal to 2). Closer to the dot, there is however no clear convergence for r\cdot T_{B}\lesssim 0.4. More precisely, the distance r_{\rm max} at which the entropy profiles reaches a maximum clearly grows when J^{\prime} goes to zero, but at some rate which is slower than T_{B}^{-1}. So, from the present data at U=2 (and V/T_{B}=17.1), there is no clear evidence of some Kondo-like spatial structure of size \sim T_{B}^{-1} in the stationary entropy profiles, as observed in the free case. It could however be that some longer times are required to achieve some profile collapse in this regime. Some further systematic investigations, as a function of U and V and J^{\prime}, would be required to elucidate this point.

## Conclusions

We have analyzed a number of numerical results concerning the steady state properties of the IRLM in the scaling regime: I-V curves, exponent b(U) and the Kondo scale T_{B}, and entropy rate \alpha. We could in particular obtain accurately the rescaled I-V curves and the entropy rate \alpha in some large range of the rescaled bias (up to V/T_{B}=100), and for several values of the interaction strength U.

We hope that this study will trigger some further investigations using analytical techniques, since our results might be compared quantitatively to field-theory results obtained thanks to the integrability of the model in the continuum limit, along the lines of what has already been done at the self-dual point Boulat, Saleur, and Schmitteckert (2008), or for the boundary sine-Gordon model Fendley, Ludwig, and Saleur (1995a, b). It would also be very interesting to elucidate qualitatively the mechanisms at work in the limit of large bias and U>0, to understand how a finite entropy rate can coexist with a vanishingly small current. Finally, our data can be used to benchmark new approximations for quantum impurity models, or when developing numerical schemes which can deal with long time evolutions in interacting problems with linear entanglement entropy growth.

## Acknowledgements

We are indebted to H. Saleur for valuable suggestions and discussions, as well as for sharing some of his unpublished notes. KB is grateful to Soltan Bidzhiev for encouragements and support in this project.

## Appendix A Details about the numerical simulations

### A.1 Initial state and Unitary evolution

The initial state is computed using a conventional DMRG procedure, and the sweeps are stopped when the energy variation becomes smaller than \Delta E=10^{-10}. During this initial state calculation the level of MPS truncation is the same as during the time evolution, and is determined by some maximum value \delta for the discarded weight, typically set to 10^{-7}.

Most of the data are calculated for chains of length between N=120 and 200. Unless mentioned otherwise, the time evolution is performed with time steps \tau=0.2 and we used the MPO-based approximation to \exp(-iH\tau) that is noted W^{\rm II} in Ref. Zaletel et al. (2015). As explained in this previous work, one can combine several W^{\rm II} to get a smaller error:

W^{\rm II}(\tau_{1})W^{\rm II}(\tau_{2})\cdots W^{\rm II}(\tau_{n})=\exp(H\tau% )+\mathcal{O}(\tau^{p}) | (10) |

The simplest case is to use n=2 steps to reduce the error to \mathcal{O}(\tau^{3}), and one solution is:

\displaystyle\tau_{1} | \displaystyle= | \displaystyle\frac{1+i}{2}\tau | |||

\displaystyle\tau_{2} | \displaystyle= | \displaystyle\frac{1-i}{2}\tau. | (11) |

We have computed two other solutions, corresponding to p=4 and p=5. These were mentioned in Ref. Zaletel et al. (2015), but not given explicitly. The first one requires n=4 and is given by

\displaystyle\tau_{1} | \displaystyle= | \displaystyle\frac{1}{4}\left(-\frac{1+i}{\sqrt{3}}+1-i\right)\tau | |||

\displaystyle\tau_{2} | \displaystyle= | \displaystyle i\tau_{1} | |||

\displaystyle\tau_{3} | \displaystyle= | \displaystyle-i\bar{\tau}_{1} | |||

\displaystyle\tau_{4} | \displaystyle= | \displaystyle\bar{\tau}_{1}. | (12) |

The next one, with an error scaling as \mathcal{O}(\tau^{p=5}) requires n=7 steps. It was obtained numerically using Mathematica:

\displaystyle\tau_{1}/\tau=0.2588533986109182+0.0447561340111419i, | ||

\displaystyle\tau_{2}/\tau=-0.0315468581488038+0.2491190542755632i, | ||

\displaystyle\tau_{3}/\tau=0.1908290521106672-0.2318537492321061i, | ||

\displaystyle\tau_{4}/\tau=0.1637288148544367, | ||

\displaystyle\tau_{5}=\bar{\tau_{3}}, | ||

\displaystyle\tau_{6}=\bar{\tau_{2}}, | ||

\displaystyle\tau_{7}=\bar{\tau_{1}}. |

It should be noted that these solutions are not unique, but we have selected the ones where the error terms, of the order \mathcal{O}(\tau^{p}), have the smallest prefactors. We have checked the precision of these three different Trotter schemes on a small free-fermion system with 6 spins, and compared the results for the current at time t=100 against the exact free-fermion solution. The results are shown in Fig. 13. As expected, the total error scales as \tau^{p-1}, due to the fact that the total number of steps is t/\tau and each step contributes an amount of order \mathcal{O}(\tau^{p}) to the error.

### A.2 Matrix truncations

During the initial DMRG sweeps (initial state calculation) and during the time evolution, the bond dimensions in the MPS are controlled using a cut off equal to \delta on the discarded weight (the sum of the discarded Schmidt values should be equal to or smaller than \delta). The effect of the truncation parameter \delta is illustrated in Figs. 14 and 15. In this case a value \delta=10^{-6} would provide a sufficient precision on the current (lower panel), but it should be noted that obtaining a precise value for the entropy rate (upper panel of Fig. 14) requires working with a smaller \delta. For smaller J^{\prime}, as shown in Fig. 15, one has to use \delta=10^{-8} in order to get some accurate estimate of the entropy rate. In that case the bond dimension can exceed 1000 at time \simeq 40.

## Appendix B Steady state in the free-fermion case

### B.1 Steady current

In the free-fermion case (U=0), the transmission and reflexion coefficients \mathcal{T}(k) and \mathcal{R}(k) for an incident fermion with momentum k have been computed by Branschadel et al. Branschädel et al. (2010):

\displaystyle\mathcal{T}(\epsilon) | \displaystyle= | \displaystyle\frac{1-\left(\epsilon/(2J)\right)^{2}}{1+\epsilon^{2}(J^{2}-2J^{% \prime 2})/\left(4J^{\prime 4}\right)} | (14) | ||

\displaystyle\mathcal{R}(\epsilon) | \displaystyle= | \displaystyle 1-\mathcal{T}(\epsilon) | (15) | ||

\displaystyle\epsilon(k) | \displaystyle= | \displaystyle-2J\cos(k) | (16) |

Combined with a Landauer approach Landauer (1957); Imry and Landauer (1999); Nazarov and Blanter (2009) this gives the steady current:

I=\int_{k_{F}^{-}}^{k_{F}^{+}}\frac{dk}{2\pi}\mathcal{T}(\epsilon(k))v(k) | (17) |

where the Fermi momenta in both leads are related to the voltage bias through \epsilon(k_{F}^{\pm})=\pm V/2 and the group velocity v(k) is, by definition, v(k)=\frac{\partial\epsilon(k)}{\partial k}. Changing the integration variable from k to \epsilon gives the standard result :

I=\int_{-V/2}^{V/2}\frac{d\epsilon}{2\pi}\mathcal{T}(\epsilon). | (18) |

For the IRLM, using Eq. 14, the integral gives:

\displaystyle 2\pi I | \displaystyle= | \displaystyle-V\frac{J^{\prime 4}}{a} | (19) | ||

\displaystyle+4J^{\prime 2}\frac{\left(1-J^{\prime 2}\right)^{2}}{a^{\frac{3}{% 2}}}\arctan\left(\frac{\sqrt{a}}{4J^{\prime 2}}V\right) | |||||

\displaystyle a | \displaystyle= | \displaystyle 1-2J^{\prime 2} |

where we assumed J=1, V<4 and J^{\prime 2}<1. This function is plotted in Fig. 16 for three different values of J^{\prime}. In the (scaling) limit where J^{\prime}\ll 1 we define T_{B}=(J^{\prime})^{2} and the Eq. B.1 simplifies to

\frac{2\pi I}{T_{B}}=4\arctan\left(\frac{V}{4T_{B}}\right), | (20) |

in agreement with Ref. Boulat, Saleur, and Schmitteckert (2008). Note that the expression I(t) of the current at finite time and U=0 can be found in the appendix of Ref. Vinkler-Aviv, Schiller, and Anders (2014). It shows damped oscillations at frequency T_{\rm osc}=4\pi/V, a relaxation time \mathcal{O}(J^{\prime 2}), and converges to Eq. 20 when t\to\infty.

### B.2 Density drop across the dot

We make the approximation that the fermions are point-like particles which propagate ballistically in leads, at some group velocity v(k). This is a semi-classical approximation (called hydrodynamical approximation in Ref. Antal, Krapivsky, and Rákos (2008)) where each particle has a well defined position and momentum. In this approximation, each lead becomes homogeneous in the steady regime and the system is described by some occupation numbers n(k)^{R/L} in both leads. Taking into account the initial momentum distributions and the scattering on the dot, we obtain:

\displaystyle n(k)^{L} | \displaystyle= | \displaystyle\left\{\begin{array}[]{cl}0,&{\rm if}\;\;k^{+}_{F}<k\\ 1,&{\rm if}\;\;-k_{F}^{-}<k<k^{+}_{F}\\ \mathcal{R}(k),&{\rm if}\;\;-k_{F}^{+}<k<-k^{-}_{F}\\ 0,&{\rm if}\;\;k<-k^{+}_{F}\\ \end{array}\right. | (21) | ||

\displaystyle n(k)^{R} | \displaystyle= | \displaystyle\left\{\begin{array}[]{cl}0,&{\rm if}\;\;k^{+}_{F}<k\\ \mathcal{T}(k),&{\rm if}\;\;k_{F}^{-}<k<k^{+}_{F}\\ 1,&{\rm if}\;\;-k_{F}^{-}<k<k^{-}_{F}\\ 0,&{\rm if}\;\;k<-k^{-}_{F}.\\ \end{array}\right. | (22) |

The total density in each lead is then obtained by integrating the distributions above. Using the fact that k_{F}^{+}+k_{F}^{-}=\pi and \mathcal{R}(k)+\mathcal{T}(k)=1 we find:

\displaystyle\rho^{L} | \displaystyle= | \displaystyle\frac{1}{2}+\int_{k_{F}^{-}}^{k_{F}^{+}}\frac{dk}{2\pi}\mathcal{R% }(\epsilon(k)) | (23) | ||

\displaystyle\rho^{R} | \displaystyle= | \displaystyle\frac{1}{2}-\int_{k_{F}^{-}}^{k_{F}^{+}}\frac{dk}{2\pi}\mathcal{R% }(\epsilon(k)). | (24) |

These densities describe the parts of the leads that are sufficiently far from the dot (|r|\gg 1), and at times that are sufficiently long (t\gg r), so that the growing quasi-steady region has reached r (and -r). On a finite system we should additionally require 2Jt\lesssim N/2 (2J being the fastest group velocity).

We thus expect some density drop

\rho^{L}-\rho^{R}=2\int_{k_{F}^{-}}^{k_{F}^{+}}\frac{dk}{2\pi}\mathcal{R}(% \epsilon(k)) | (25) |

across the dot.

### B.3 Entropy rate

In analogy with the Landauer approach for the current, the stationary entropy rate \alpha can be obtained analytically. Since each incident particle at energy \epsilon contributes to the entropy by an amount \delta S=-\mathcal{T}(\epsilon)\ln\mathcal{T}(\epsilon)-\mathcal{R}(\epsilon)% \ln\mathcal{R}(\epsilon) (its wave packet is split into a reflected part and a transmitted part), we get:

\alpha=-\frac{1}{2\pi}\int_{-V/2}^{V/2}d\epsilon\left[\mathcal{T}(\epsilon)\ln% \mathcal{T}(\epsilon)+\mathcal{R}(\epsilon)\ln\mathcal{R}(\epsilon)\right]. | (26) |

The result above has been checked numerically against numerical solution of dynamics for the free-fermion problem.

In the scaling regime J^{\prime} is small, the energy \epsilon are small (because the bias V is small), but the ratio \epsilon/J^{\prime 2} is of order one. In this limit the transmission coefficient (Eq. 14) becomes:

\mathcal{T}(x)=\frac{1}{1+x^{2}} | (27) |

where x=\epsilon/(2T_{B}) and T_{B}=J^{\prime 2}. With a change of variable the entropy rate (Eq. 26) can be expressed as:

\displaystyle\alpha(v) | \displaystyle= | \displaystyle\frac{2T_{B}}{2\pi}\int_{-v/4}^{v/4}dx\left(\frac{1}{1+x^{2}}\ln(% 1+x^{2})\right. | (28) | ||

\displaystyle\left.+\frac{x^{2}}{1+x^{2}}\ln\left(\frac{1+x^{2}}{x^{2}}\right)% \right), |

where v=V/T_{B}. The integral can be computed explicitly and the final result is

\displaystyle\frac{\alpha(v)}{T_{B}} | \displaystyle= | \displaystyle\frac{2}{\pi}\left[2\left(1+\ln(v/4)\right)\arctan\left(v/4\right% )\right. | (29) | ||

\displaystyle+\frac{1}{4}\,v\ln\left(v^{2}+16\right)-\frac{1}{2}v\ln(v) | |||||

\displaystyle\left.-i{\rm Li}_{2}\left(-\frac{iv}{4}\right)+i{\rm Li}_{2}\left% (\frac{iv}{4}\right)\right] |

where {\rm Li}_{2} is the the polylogarithm of index 2. This quantity tends to a constant at large bias:

\alpha(v\to\infty)/T_{B}=2. | (30) |

And at low bias we have:

\frac{\pi}{T_{B}}\alpha(v\to 0)=\left(\frac{5}{288}-\frac{1}{48}\ln(v/4)\right% )v^{3}+\mathcal{O}\left(v^{5}\right). | (31) |

## References

- Polkovnikov et al. (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
- Eisert, Friesdorf, and Gogolin (2011) J. Eisert, M. Friesdorf, and C. Gogolin, Nat. Phys. 11, 124 (2011).
- Schollwöck (2011) U. Schollwöck, Ann. Phys.ics 326, 96 (2011).
- Hewson (1993) A. C. Hewson, The Kondo Problem to Heavy Fermions (Cambridge University Press, 1993).
- Wiegmann and Finkel’shtein (1978) P. B. Wiegmann and A. M. Finkel’shtein, Sov. Phys. JETP 48, 102 (1978).
- Filyov and Wiegmann (1980) V. M. Filyov and P. B. Wiegmann, Phys. Lett. A 76, 283 (1980).
- Boulat, Saleur, and Schmitteckert (2008) E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. Lett. 101, 140601 (2008).
- Fendley, Ludwig, and Saleur (1995a) P. Fendley, A. W. W. Ludwig, and H. Saleur, Phys. Rev. Lett. 74, 3005 (1995a).
- Fendley, Ludwig, and Saleur (1995b) P. Fendley, A. W. W. Ludwig, and H. Saleur, Phys. Rev. B 52, 8934 (1995b).
- Mehta and Andrei (2006) P. Mehta and N. Andrei, Phys. Rev. Lett. 96, 216802 (2006).
- Gobert et al. (2005) D. Gobert, C. Kollath, U. Schollwöck, and G. Schütz, Phys. Rev. E 71, 036102 (2005).
- Sabetta and Misguich (2013) T. Sabetta and G. Misguich, Phys. Rev. B 88, 245114 (2013).
- Schmitteckert (2004) P. Schmitteckert, Phys. Rev. B 70, 121302 (2004).
- Schneider and Schmitteckert (2006) G. Schneider and P. Schmitteckert, arXiv:cond-mat/0601389 (2006).
- Vidal (2004) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
- White and Feiguin (2004) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
- Daley et al. (2004) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, J. Stat. Mech. 2004, P04005 (2004).
- Laflorencie (2016) N. Laflorencie, Phys. Rep. 646, 1 (2016).
- ite (n 20) ITensor Library, http://itensor.org (version 2.0).
- Zaletel et al. (2015) M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore, and F. Pollmann, Phys. Rev. B 91, 165112 (2015).
- Eisler and Peschel (2012) V. Eisler and I. Peschel, EPL 99, 20001 (2012).
- Kennes, Meden, and Vasseur (2014) D. M. Kennes, V. Meden, and R. Vasseur, Phys. Rev. B 90, 115101 (2014).
- Vasseur and Saleur (2017) R. Vasseur and H. Saleur, SciPost Physics 3, 001 (2017).
- Calabrese and Cardy (2007) P. Calabrese and J. Cardy, J. Stat. Mech. 2007, P10004 (2007).
- Stéphan and Dubail (2011) J.-M. Stéphan and J. Dubail, J. Stat. Mech. 2011, P08019 (2011).
- Klich and Levitov (2009) I. Klich and L. Levitov, Phys. Rev. Lett. 102, 100502 (2009).
- Song et al. (2012) H. F. Song, S. Rachel, C. Flindt, I. Klich, N. Laflorencie, and K. Le Hur, Phys. Rev. B 85, 035409 (2012).
- Song et al. (2011) H. F. Song, C. Flindt, S. Rachel, I. Klich, and K. Le Hur, Phys. Rev. B 83, 161408 (2011).
- Karrasch et al. (2010) C. Karrasch, M. Pletyukhov, L. Borda, and V. Meden, Phys. Rev. B 81, 125122 (2010).
- Vinkler-Aviv, Schiller, and Anders (2014) Y. Vinkler-Aviv, A. Schiller, and F. B. Anders, Phys. Rev. B 90, 155110 (2014).
- Freton, Boulat, and Saleur (2013) L. Freton, E. Boulat, and H. Saleur, Nucl. Phys. B 874, 279 (2013).
- Branschädel et al. (2010) A. Branschädel, E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. B 82, 205414 (2010).
- Landauer (1957) R. Landauer, IBM J. Res. Dev. 1, 223 (1957).
- Imry and Landauer (1999) Y. Imry and R. Landauer, Rev. Mod. Phys. 71, S306 (1999).
- Nazarov and Blanter (2009) Y. V. Nazarov and Y. Blanter, Quantum transport - Introduction to Nanoscience (Cambridge University Press, 2009).
- Antal, Krapivsky, and Rákos (2008) T. Antal, P. L. Krapivsky, and A. Rákos, Phys. Rev. E 78, 061115 (2008).