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

Kemal Bidzhiev, Grégoire Misguich

<figure><img src="image_1.png" /></figure>

## To cite this version:

Kemal Bidzhiev, Grégoire Misguich. Out-of-equilibrium dynamics in a quantum impurity model: Numerics for particle transport and entanglement entropy. Physical Review B, 2017, 96 (19), pp.195117. 10.1103/PhysRevB.96.195117. hal-01645813

HAL Id: hal-01645813

https://hal.science/hal-01645813

Submitted on 18 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Out-of-equilibrium dynamics in a quantum impurity model: Numerics for particle transport and entanglement entropy

Kemal Bidzhiev, Grégoire Misguich

► To cite this version:

Kemal Bidzhiev, Grégoire Misguich. Out-of-equilibrium dynamics in a quantum impurity model: Numerics for particle transport and entanglement entropy. Physical Review B, 2017, 96 (19), pp.195117. 10.1103/PhysRevB.96.195117. hal-01645813

HAL Id: hal-01645813

https://hal.science/hal-01645813v1

Submitted on 18 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

Kemal Bidzhiev and Grégoire Misguich

Institut de Physique Théorique, Université Paris Saclay, CEA, CNRS, F-91191 Gif-sur-Yvette, France

(Received 19 July 2017; revised manuscript received 27 October 2017; published 8 November 2017)

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.

DOI: 10.1103/PhysRevB.96.195117

## I. INTRODUCTION

The dynamical properties of out-of-equilibrium quantum many-body systems are a major topic in condensed-matter physics [1,2]. 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 [3], 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 [4] 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) [5]. 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 [6], several remarkable exact results have been obtained concerning nonequilibrium steady states, such as the current-voltage (I-V) characteristic for some special (so-called “self-dual”) value of the interaction strength [7] (see also [8–10]). 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 quasisteady 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-1/2 chains [11,12], or impurity models [7,13,14]. These simulations have been performed using the time-evolving block decimation and time-dependent density matrix renormalization group (DMRG) [15–17], 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 [7], 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 [18]). 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,
<figure><img src="image_2.png" /></figure>

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

Kemal Bidzhiev and Grégoire Misguich

Institut de Physique Théorique, Université Paris Saclay, CEA, CNRS, F-91191 Gif-sur-Yvette, France

(Received 19 July 2017; revised manuscript received 27 October 2017; published 8 November 2017)

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.

DOI: 10.1103/PhysRevB.96.195117

## I. INTRODUCTION

The dynamical properties of out-of-equilibrium quantum many-body systems are a major topic in condensed-matter physics [1,2]. 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 [3], 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 [4] 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) [5]. 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 [6], several remarkable exact results have been obtained concerning nonequilibrium steady states, such as the current-voltage (I-V) characteristic for some special (so-called “self-dual”) value of the interaction strength [7] (see also [8–10]). 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 quasisteady 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-1/2 chains [11,12], or impurity models [7,13,14]. These simulations have been performed using the time-evolving block decimation and time-dependent density matrix renormalization group (DMRG) [15–17], 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 [7], 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 [18]). 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,
<figure><img src="image_4.png" /><figcaption>FIG. 1. Schematics of the interacting resonant level model. The system is prepared at $t = 0$ in the ground state of the model with an additional chemical potential equal to $+V/2$ in the left lead and equal to $-V/2$ in the right lead. For $t > 0$, the system then evolves with the bias switched to zero. In all cases, the chemical potential is zero on the dot, hence the name “resonant”. Note that, by symmetry, the average fermion occupation number $\langle c_0^\dagger c_0 \rangle$ on the dot is equal to $1/2$.</figcaption></figure>

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 B), following Ref. [7], 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 C, we present a similar analysis for the scaling of the entropy rate, and we show that it can also be analyzed in terms 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 D.

Finally, Appendix A provides some technical details concerning the numerical simulations, and Appendix B presents a few exact results in the free-fermion case (current, density, and entropy rate).

## II. MODEL AND TIME EVOLUTION

### A. Hamiltonian

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

$$H_{\text{IRLM}} = H_A + H_B + H_d, \quad (1)$$

$$H_A = -J \sum_{r=-N/2}^{N/2-1} (c_r^\dagger c_{r+1} + \text{H.c.}), \quad (2)$$

$$H_B = -J \sum_{r=1}^{N/2-1} (c_r^\dagger c_{r+1} + \text{H.c.}), \quad (3)$$

$$H_d = -J' \sum_{r=-1}^{0} (c_r^\dagger c_{r+1} + \text{H.c.}) \\ + U \sum_{r=\pm 1} \left(c_r^\dagger c_r - \frac{1}{2}\right) \left(c_0^\dagger c_0 - \frac{1}{2}\right), \quad (4)$$

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'$ (or tunneling strength) between the leads

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 2. Magnetization profile $S^z(r,t)$ at times $t = 0, 6, 12, 18, 24$. Parameters of the model: $J' = 0.3$, $V = 2.0$, $N = 254$. Top panel: $U = 2.0$, bottom: $U = 0$. In both cases, the initial magnetization $m_0$ in the bulk of the left lead is indicated by a (red) horizontal line. The orange horizontal line marks the exact bulk stationary magnetization $m$ for $U = 0$ [see Eq. (6)]. The data in the upper panel show some deviations from this noninteracting value. Note that some Friedel-type oscillations develop at long times in the vicinity of the dot, and these are much stronger in the noninteracting case.</figcaption></figure>

and the dot ($r = 0$) is much smaller than the bandwidth $4J$ of the kinetic energy in the leads, i.e., $J' \ll J$. From now we take $J = 1 = \hbar$, thus defining the unit of time and energy.

### B. Initial state

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

$$H_{\text{bias}} = \frac{1}{2} V \sum_{r=-N/2}^{N/2-1} \tanh(r/w) c_r^\dagger c_r. \quad (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_r^\dagger 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^\pm (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 [7], 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
<figure><img src="image_3.png" /><figcaption>FIG. 1. Schematics of the interacting resonant level model. The system is prepared at $t = 0$ in the ground state of the model with an additional chemical potential equal to $+V/2$ in the left lead and equal to $-V/2$ in the right lead. For $t > 0$, the system then evolves with the bias switched to zero. In all cases, the chemical potential is zero on the dot, hence the name “resonant”. Note that, by symmetry, the average fermion occupation number $\langle c_0^\dagger c_0 \rangle$ on the dot is equal to $1/2$.</figcaption></figure>

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 B), following Ref. [7], 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 C, we present a similar analysis for the scaling of the entropy rate, and we show that it can also be analyzed in terms 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 D.

Finally, Appendix A provides some technical details concerning the numerical simulations, and Appendix B presents a few exact results in the free-fermion case (current, density, and entropy rate).

## II. MODEL AND TIME EVOLUTION

### A. Hamiltonian

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

$$H_{\text{IRLM}} = H_A + H_B + H_d, \quad (1)$$

$$H_A = -J \sum_{r=-N/2}^{N/2-1} (c_r^\dagger c_{r+1} + \text{H.c.}), \quad (2)$$

$$H_B = -J \sum_{r=1}^{N/2-1} (c_r^\dagger c_{r+1} + \text{H.c.}), \quad (3)$$

$$H_d = -J' \sum_{r=-1}^{0} (c_r^\dagger c_{r+1} + \text{H.c.}) + U \sum_{r=\pm 1} \left(c_r^\dagger c_r - \frac{1}{2}\right) \left(c_0^\dagger c_0 - \frac{1}{2}\right), \quad (4)$$

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'$ (or tunneling strength) between the leads

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 2. Magnetization profile $S^z(r,t)$ at times $t = 0, 6, 12, 18, 24$. Parameters of the model: $J' = 0.3$, $V = 2.0$, $N = 254$. Top panel: $U = 2.0$, bottom: $U = 0$. In both cases, the initial magnetization $m_0$ in the bulk of the left lead is indicated by a (red) horizontal line. The orange horizontal line marks the exact bulk stationary magnetization $m$ for $U = 0$ [see Eq. (6)]. The data in the upper panel show some deviations from this noninteracting value. Note that some Friedel-type oscillations develop at long times in the vicinity of the dot, and these are much stronger in the noninteracting case.</figcaption></figure>

and the dot ($r = 0$) is much smaller than the bandwidth $4J$ of the kinetic energy in the leads, i.e., $J' \ll J$. From now we take $J = 1 = \hbar$, thus defining the unit of time and energy.

### B. Initial state

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

$$H_{\text{bias}} = \frac{1}{2} V \sum_{r=-N/2}^{N/2-1} \tanh(r/w) c_r^\dagger c_r. \quad (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_r^\dagger 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^\pm (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 [7], 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' = 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$.

### C. Unitary evolution at $t > 0$

For $t > 0$ the wave function evolves using $H_{IRLM}$, with $0 < J' < 1$, and without the voltage bias term. The calculations are performed using a time-dependent DMRG algorithm, implemented using the C++ iTensor library [19]. The evolution operator $U = \exp(-i\tau H)$ for a time step $\tau$ is approximated by a matrix-product operator (MPO) [20], using a fourth-order Trotter scheme [Eqs. (A1) and (A4)] 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 \approx 100$. Additional details about the simulations are given in Appendix A.

In the following, we focus on three quantities: the particle density (Sec. IID), the particle current (Sec. IIE), and the entanglement entropy (Sec. IIF).

### D. Particle density

The particle density is defined by $\rho(r,t) = \langle\psi(t)|c_r^\dagger 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 regions form in the center. When the time is large enough, two regions with quasihomogeneous 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)), \quad (6)$$

where $\mathcal{R}(\epsilon)$ is the reflexion coefficient for an incident particle with energy $\epsilon$ (more details in Appendix B 2). This exact value $m$ is in agreement with the magnetization that is measured numerically in the stationary region for $U = 0$ (bottom panel of Fig. 2), but slightly different from that observed in the interacting case (top panel of Fig. 2), as expected.

### E. Particle current

On a given bond, the expectation value of the current operator is $I(r,t) = 2J_r \mathrm{Im}\langle\psi(t)|c_r^\dagger 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'$ 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_r^\dagger 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}[I(-1,t) + I(0,t)]$.

<figure><img src="image_7.png" /><figcaption>FIG. 3. Top panels: evolution of the entropy $S(r,t)$ profile for $U = 0$ and 2. Bottom panels: evolution of the current profile $I(r,t)$ for $U = 0$ and 2.0. Parameters of the model: $J' = 0.3$, $V = 2.0$, $N = 254$, as in Fig. 2.</figcaption></figure>

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
the convergence to a steady state. We typically use $w = 10$. For the same reason, the initial state is prepared with $J' = 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$.

### C. Unitary evolution at $t > 0$

For $t > 0$ the wave function evolves using $H_{IRLM}$, with $0 < J' < 1$, and without the voltage bias term. The calculations are performed using a time-dependent DMRG algorithm, implemented using the C++ iTensor library [19]. The evolution operator $U = \exp(-i\tau H)$ for a time step $\tau$ is approximated by a matrix-product operator (MPO) [20], using a fourth-order Trotter scheme [Eqs. (A1) and (A4)] 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.

In the following, we focus on three quantities: the particle density (Sec. II D), the particle current (Sec. II E), and the entanglement entropy (Sec. II F).

### D. Particle density

The particle density is defined by $\rho(r,t) = \langle\psi(t)|c_r^\dagger 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 regions form in the center. When the time is large enough, two regions with quasihomogeneous 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)), \quad (6)$$

where $\mathcal{R}(\epsilon)$ is the reflexion coefficient for an incident particle with energy $\epsilon$ (more details in Appendix B 2). This exact value $m$ is in agreement with the magnetization that is measured numerically in the stationary region for $U = 0$ (bottom panel of Fig. 2), but slightly different from that observed in the interacting case (top panel of Fig. 2), as expected.

### E. Particle current

On a given bond, the expectation value of the current operator is $I(r,t) = 2J_r \text{Im}\langle\psi(t)|c_r^\dagger 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'$ 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_r^\dagger 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}[I(-1,t) + I(0,t)]$.

<figure><img src="image_6.png" /><figcaption>FIG. 3. Top panels: evolution of the entropy $S(r,t)$ profile for $U = 0$ and 2. Bottom panels: evolution of the current profile $I(r,t)$ for $U = 0$ and 2.0. Parameters of the model: $J' = 0.3$, $V = 2.0$, $N = 254$, as in Fig. 2.</figcaption></figure>

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
<figure><img src="image_8.png" /><figcaption>FIG. 4. Top panel: current $I(t)$. Bottom panel: entanglement entropy $S(t)$ between the left lead and the rest of the system. The horizontal line in the top panel is the exact result at the self-dual point (see text), derived in Ref. [7]. Parameters: $U = 2.0$, $J' = 0.3$, $V = 2.0$, $N = 254$.</figcaption></figure>

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_{\text{osc}}$. It was argued in Ref. [14] that this period is simply given by the bias: $T_{\text{osc}} = 4\pi/V$. Our data are in agreement with this result, from $U=0$ to large values of $U$.

### F. 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. [21]. 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 [22,23]. A logarithmic growth is in fact generic for local quenches

in critical one-dimensional systems [21,24,25]. 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 nonzero bias) [11,12].

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' = 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 [26–28]. 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 C 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$.<fnref n="1" />

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. (B13). Although this classical picture with noninteracting particles does not apply to the interacting case, the numerical simulations show

<footnote n="1">¹The 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. (B7). Note that the actual distribution is in fact Poissonian.</footnote>
<figure><img src="image_7.png" /><figcaption>FIG. 4. Top panel: current $I(t)$. Bottom panel: entanglement entropy $S(t)$ between the left lead and the rest of the system. The horizontal line in the top panel is the exact result at the self-dual point (see text), derived in Ref. [7]. Parameters: $U = 2.0$, $J' = 0.3$, $V = 2.0$, $N = 254$.</figcaption></figure>

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_{\text{osc}}$. It was argued in Ref. [14] that this period is simply given by the bias: $T_{\text{osc}} = 4\pi/V$. Our data are in agreement with this result, from $U=0$ to large values of $U$.

### F. 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. [21]. 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 [22,23]. A logarithmic growth is in fact generic for local quenches

in critical one-dimensional systems [21,24,25]. 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 nonzero bias) [11,12].

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' = 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 [26–28]. 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 C 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$

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. (B13). Although this classical picture with noninteracting particles does not apply to the interacting case, the numerical simulations show

<footnote>$^1$The 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. (B7). Note that the actual distribution is in fact Poissonian.</footnote>
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. [7–9], 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 D.

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.

#### A. 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'$, $V$, and $U$.$^2$

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:

$$H_A^c = i v_F \int_{-\infty}^{0} dr [\psi_L^\dagger(r) \partial_r \psi_L(r) - \psi_R^\dagger(r) \partial_r \psi_R(r)], \quad (7)$$

$$H_B^c = i v_F \int_{0}^{\infty} dr [\psi_L^\dagger(r) \partial_r \psi_L(r) - \psi_R^\dagger(r) \partial_r \psi_R(r)]. \quad (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

$$H_d^c = -J'_c(\psi_L(0) + \psi_R(0))d^\dagger + \text{H.c.} \\ + U_c(\psi_L(0)^\dagger\psi_L(0) + \psi_R(0)^\dagger\psi_R(0)) \\ \times (d^\dagger d - \frac{1}{2}), \quad (9)$$

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)$ [7]. 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$ [22,23].

In Sec. II F, 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'^{1/(1-\alpha)}$, 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} + (\frac{1}{2} - \frac{U_c}{2\pi})^2$, where $U_c$ is the interaction strength in the continuum limit [7]. 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 infrared fixed point, but there is no known method to do this exactly. The analytical result of Ref. [7] 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 [7]. 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$ [7]. This behavior will be checked numerically in detail in the next section (Sec. III B). 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'$.

#### B. Steady current

Figure 5 shows the stationary current as a function of the bias, for a few values of $U$ and $J' = 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'$. 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$ Then, the scale $T_B$ is defined

<footnote>$^2$In practice we take $J = 1$ and restrict our numerical simulations to $V \lesssim 2$, $J' \lesssim 0.3$, and $U \lesssim 2.5$. Finite-$J'$ corrections start to be visible for $J' = 0.5$. As for $U$, our data at $U = 4$ are clearly outside the scaling regime.</footnote>

<footnote>$^3$The value $J' = 0.08$ appears to offer a good compromise to estimate $b$ in our simulations. Indeed, we need a small $J'$ to be in the scaling regime, but the time to reach the steady regime (which increases when $J'$ decreases) should also not become too large compared to $L/v_F$.</footnote>
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. [7–9], 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 D.

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.

#### A. 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'$, $V$, and $U$.$^2$

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:

$$H_A^c = i v_F \int_{-\infty}^{0} dr [\psi_L^\dagger(r) \partial_r \psi_L(r) - \psi_R^\dagger(r) \partial_r \psi_R(r)], \quad (7)$$

$$H_B^c = i v_F \int_0^\infty dr [\psi_L^\dagger(r) \partial_r \psi_L(r) - \psi_R^\dagger(r) \partial_r \psi_R(r)]. \quad (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

$$H_d^c = -J_c'(\psi_L(0) + \psi_R(0))d^\dagger + \text{H.c.} \\ + U_c(\psi_L(0)^\dagger\psi_L(0) + \psi_R(0)^\dagger\psi_R(0)) \\ \times (d^\dagger d - \frac{1}{2}), \quad (9)$$

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)$ [7]. 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$ [22,23].

In Sec. II F, 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'^{1/(1-\alpha)}$, 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} + (\frac{1}{2} - \frac{U_c}{2\pi})^2$, where $U_c$ is the interaction strength in the continuum limit [7]. 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 infrared fixed point, but there is no known method to do this exactly. The analytical result of Ref. [7] 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 [7]. 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$ [7]. This behavior will be checked numerically in detail in the next section (Sec. III B). 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'$.

#### B. Steady current

Figure 5 shows the stationary current as a function of the bias, for a few values of $U$ and $J' = 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'$. 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$ Then, the scale $T_B$ is defined

<footnote>$^2$In practice we take $J = 1$ and restrict our numerical simulations to $V \lesssim 2$, $J' \lesssim 0.3$, and $U \lesssim 2.5$. Finite-$J'$ corrections start to be visible for $J' = 0.5$. As for $U$, our data at $U = 4$ are clearly outside the scaling regime.</footnote>

<footnote>$^3$The value $J' = 0.08$ appears to offer a good compromise to estimate $b$ in our simulations. Indeed, we need a small $J'$ to be in the scaling regime, but the time to reach the steady regime (which increases when $J'$ decreases) should also not become too large compared to $L/v_F$.</footnote>
<figure><img src="image_9.png" /><figcaption>FIG. 5. Stationary current *I* as a function of *V* for *J*' = 0.08, for a few selected values of *U*. A log-log scale is used to show the power-law behavior of the current at large *V*/*T*<sub>B</sub> (*I* ~ *V*<sup>-b</sup>), and to extract the associated exponent *b*(*U*). The values *b*(*U*) extracted by these fits are displayed in Fig. 6.</figcaption></figure>

as $T_B(J',U) = (J')^{\frac{1}{1-b(U)}} = (J')^{\frac{2}{1+b(U)}}$.<fnref n="4" /> From the analysis of the model in the continuum, we know that the exponent *b* should reach a maximum value *b* = 1/2 (equivalent to *h* = 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'$ (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.* [7], but thanks to longer simulations and larger systems, the present

<figure><img src="image_10.png" /><figcaption>FIG. 6. Exponent $b(U)$ as a function of the interaction strength $U$, obtained by fitting the steady current $I$ to $I \sim V^{-b}$, at large $V/T_B$ (see Fig. 5). The full (black) line is $b = 2U/\pi$, the result of some small-$U$ expansion [29,30].</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 7. Rescaled current $I/T_B$ as a function of the rescaled bias $V/T_B$, for different values of $U$. The colors label the values of $U$, while the symbol shapes encode $J'$ (see the legend of Fig. 8 for details). For a given $U$, the results obtained for different values of $J'$ approximately collapse onto a single curve, as expected in the scaling regime. The red line is the theoretical result for the self-dual point [7]. The numerical data for $U = 2.0$ appear to be in very good agreement with this theoretical curve. The black line is the exact result for $U = 0$ [Eq. (B7)]. $T_B$ is defined as $T_B = (J')^{\frac{1}{1-b(U)}}$ where the exponent $h(U)$ is determined from the behavior of $I$ at large $V/T_B$ (see text and Fig. 6).</figcaption></figure>

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 \le 0$ the current is monotonically increasing), and has already been discussed in [7]. For small values of $U$, the exponent $b$ describing the current suppression at large bias has been computed using some functional renormalization group method [29] or with some self-consistent conserving approximation [30]. Using our convention for the lattice model, their result reads as $b = 2U/\pi + 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'$ (0.08) and finite voltage $V$ (up to $\approx 2$).

### C. Entropy rate

We estimate the steady entropy rate $\alpha$, defined in Sec. II F, by fitting the long-time part of the entanglement entropy data $S(t)$.<fnref n="5" /> 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

<footnote n="4">⁴We adopt the convention where the numerical prefactor in the definition of $T_B$ is set to 1.</footnote>

<footnote n="5">⁵Like for the current $I(t)$, the entropy $S(t)$ often shows some oscillatory part, ~$\cos(Vt/2 + \text{cst})$ (see Fig. 15), and we include such a term in the fitting function in order to improve the precision on $\alpha$.</footnote>
<figure><img src="image_8.png" /><figcaption>FIG. 5. Stationary current *I* as a function of *V* for *J*' = 0.08, for a few selected values of *U*. A log-log scale is used to show the power-law behavior of the current at large *V*/*T*<sub>B</sub> (*I* ~ *V*<sup>-b</sup>), and to extract the associated exponent *b*(*U*). The values *b*(*U*) extracted by these fits are displayed in Fig. 6.</figcaption></figure>

as $T_B(J',U) = (J')^{\frac{1}{1-b(U)}} = (J')^{\frac{2}{1+b(U)}}$.<fnref n="4" /> From the analysis of the model in the continuum, we know that the exponent *b* should reach a maximum value *b* = 1/2 (equivalent to *h* = 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'$ (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.* [7], but thanks to longer simulations and larger systems, the present

<figure><img src="image_9.png" /><figcaption>FIG. 6. Exponent $b(U)$ as a function of the interaction strength $U$, obtained by fitting the steady current $I$ to $I \sim V^{-b}$, at large $V/T_B$ (see Fig. 5). The full (black) line is $b = 2U/\pi$, the result of some small-$U$ expansion [29,30].</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 7. Rescaled current $I/T_B$ as a function of the rescaled bias $V/T_B$, for different values of $U$. The colors label the values of $U$, while the symbol shapes encode $J'$ (see the legend of Fig. 8 for details). For a given $U$, the results obtained for different values of $J'$ approximately collapse onto a single curve, as expected in the scaling regime. The red line is the theoretical result for the self-dual point [7]. The numerical data for $U = 2.0$ appear to be in very good agreement with this theoretical curve. The black line is the exact result for $U = 0$ [Eq. (B7)]. $T_B$ is defined as $T_B = (J')^{\frac{1}{1-b(U)}}$ where the exponent $b$ is determined from the behavior of $I$ at large $V/T_B$ (see text and Fig. 6).</figcaption></figure>

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 \le 0$ the current is monotonically increasing), and has already been discussed in [7]. For small values of $U$, the exponent $b$ describing the current suppression at large bias has been computed using some functional renormalization group method [29] or with some self-consistent conserving approximation [30]. Using our convention for the lattice model, their result reads as $b = 2U/\pi + 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'$ (0.08) and finite voltage $V$ (up to $\simeq 2$).

### C. Entropy rate

We estimate the steady entropy rate $\alpha$, defined in Sec. II F, by fitting the long-time part of the entanglement entropy data $S(t)$.<fnref n="5" /> 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

<footnote n="4">⁴We adopt the convention where the numerical prefactor in the definition of $T_B$ is set to 1.</footnote>

<footnote n="5">⁵Like for the current $I(t)$, the entropy $S(t)$ often shows some oscillatory part, ~$\cos(Vt/2 + \text{cst})$ (see Fig. 15), and we include such a term in the fitting function in order to improve the precision on $\alpha$.</footnote>
<figure><img src="image_12.png" /><figcaption>FIG. 8. Same as Fig. 7 (rescaled current $I/T_B$ as a function of the rescaled bias $V/T_B$), with a zoom on the low-bias region. The colors label the values of $U$, while the symbol shapes encode $J'$ (see the legend).</figcaption></figure>

given value $U$ but for different values of $J'$ 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' = 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.

<figure><img src="image_13.png" /><figcaption>FIG. 9. Rescaled entropy rate, $\alpha/T_B$, as a function of $V/T_B$, for several values of $U$. The colors label the values of $U$, while the symbol shapes encode $J'$ (see the legend of Fig. 8 for details). The black line is the exact free-fermion result [Eq. (B16)].</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 10. Same as Fig. 9, with a zoom on the low-bias region. Top panel: $U = -0.1, 0, \text{ and } 3$. Bottom panel: $U = 0.2, 1, \text{ and } 2$. The colors encode the values of $U$, while the symbol shapes label the values of $J'$, as indicated in the legend of Fig. 8. The black line (top panel) is the exact free-fermion result in the limit $J' \to 0$ [Eq. (B16)].</figcaption></figure>

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

### D. 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. [23,31] concerning the IRLM, and Ref. [18] for a more general review. While these studies investigated the entanglement in the ground state
<figure><img src="image_11.png" /><figcaption>FIG. 8. Same as Fig. 7 (rescaled current $I/T_B$ as a function of the rescaled bias $V/T_B$), with a zoom on the low-bias region. The colors label the values of $U$, while the symbol shapes encode $J'$ (see the legend).</figcaption></figure>

given value $U$ but for different values of $J'$ 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' = 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.

<figure><img src="image_12.png" /><figcaption>FIG. 9. Rescaled entropy rate, $\alpha/T_B$, as a function of $V/T_B$, for several values of $U$. The colors label the values of $U$, while the symbol shapes encode $J'$ (see the legend of Fig. 8 for details). The black line is the exact free-fermion result [Eq. (B16)].</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 10. Same as Fig. 9, with a zoom on the low-bias region. Top panel: $U = -0.1, 0, \text{ and } 3$. Bottom panel: $U = 0.2, 1, \text{ and } 2$. The colors encode the values of $U$, while the symbol shapes label the values of $J'$, as indicated in the legend of Fig. 8. The black line (top panel) is the exact free-fermion result in the limit $J' \to 0$ [Eq. (B16)].</figcaption></figure>

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

### D. 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. [23,31] concerning the IRLM, and Ref. [18] for a more general review. While these studies investigated the entanglement in the ground state
<figure><img src="image_15.png" /></figure>

FIG. 11. Top: entanglement entropy profiles for $U = 0$, plotted as a function of the rescaled distance $r/\xi$, with $\xi = T_B^{-1}$ and $T_B = J^2$. Some constant value has been subtracted from each profile, to allow the data to (approximately) fall onto a single time-independent curve. We only show here the right part of the profile ($r > 0$), the other side ($r < 0$) being symmetric. Bottom: same data plotted as a function of the bare distance $r$. System size is $N = 6000$ and time is $t = 1200$, values for $V(J')$: $V(0.25) = 0.5$, $V(0.1) = 0.08$, $V(0.075) = 0.045$, $V(0.05) = 0.02$ with ratio $T_B^{-1}(0.05):T_B^{-1}(0.07):T_B^{-1}(0.1):T_B^{-1}(0.25) = 25:4:2.25:1$. In all cases, the rescaled bias is $V/T_B = 8.0$.

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 12. As discussed in Sec. II F, the global shape of the profile is approximately triangular, with some maximum value that grows as $\alpha t$ with a rate $\alpha \sim 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

<figure><img src="image_16.png" /></figure>

FIG. 12. Same as Fig. 11, but for $U = 2.0$. System size is $N = 402$ and $t = 60$. Values for $V(J')$: $V(0.12) = 1.0$, $V(0.07) = 0.5$, $V(0.04) = 0.25$, $V(0.02) = 0.1$ with ratio $T_B^{-1}(0.02):T_B^{-1}(0.04):T_B^{-1}(0.07):T_B^{-1}(1.2) = 10:4:2:1$. In all cases, the rescaled bias is $V/T_B = 17.1$.

values of $J'$ (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'$ (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_{max}$ at which the entropy profiles reach a maximum clearly grows when $J'$ 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-type 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'$, would be required to elucidate this point.

## IV. 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 [7], or for the boundary sine-Gordon model [8,9]. 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.

## ACKNOWLEDGMENTS

We are indebted to H. Saleur for valuable suggestions and discussions, as well as for sharing some of his unpublished notes.
<figure><img src="image_14.png" /></figure>

FIG. 11. Top: entanglement entropy profiles for $U = 0$, plotted as a function of the rescaled distance $r/\xi$, with $\xi = T_B^{-1}$ and $T_B = J^2$. Some constant value has been subtracted from each profile, to allow the data to (approximately) fall onto a single time-independent curve. We only show here the right part of the profile ($r > 0$), the other side ($r < 0$) being symmetric. Bottom: same data plotted as a function of the bare distance $r$. System size is $N = 6000$ and time is $t = 1200$, values for $V(J')$: $V(0.25) = 0.5$, $V(0.1) = 0.08$, $V(0.075) = 0.045$, $V(0.05) = 0.02$ with ratio $T_B^{-1}(0.05): T_B^{-1}(0.07): T_B^{-1}(0.1): T_B^{-1}(0.25) = 25:4:2.25:1$. In all cases, the rescaled bias is $V/T_B = 8.0$.

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 12. As discussed in Sec. II F, the global shape of the profile is approximately triangular, with some maximum value that grows as $\alpha t$ with a rate $\alpha \sim 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

<figure><img src="image_15.png" /></figure>

FIG. 12. Same as Fig. 11, but for $U = 2.0$. System size is $N = 402$ and $t = 60$. Values for $V(J')$: $V(0.12) = 1.0$, $V(0.07) = 0.5$, $V(0.04) = 0.25$, $V(0.02) = 0.1$ with ratio $T_B^{-1}(0.02): T_B^{-1}(0.04): T_B^{-1}(0.07): T_B^{-1}(1.2) = 10:4:2:1$. In all cases, the rescaled bias is $V/T_B = 17.1$.

values of $J'$ (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'$ (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_{max}$ at which the entropy profiles reach a maximum clearly grows when $J'$ 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-type 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'$, would be required to elucidate this point.

## IV. 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 [7], or for the boundary sine-Gordon model [8,9]. 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.

## ACKNOWLEDGMENTS

We are indebted to H. Saleur for valuable suggestions and discussions, as well as for sharing some of his unpublished notes.
# APPENDIX A: DETAILS ABOUT THE NUMERICAL SIMULATIONS

## 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^{\text{II}}$ in Ref. [20]. As explained in this previous work, one can combine several $W^{\text{II}}$ to get a smaller error:

$$W^{\text{II}}(\tau_1)W^{\text{II}}(\tau_2)\dots W^{\text{II}}(\tau_n) = \exp(H\tau) + O(\tau^p). \quad (A1)$$

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

$$\tau_1 = \frac{1+i}{2}\tau, \quad \tau_2 = \frac{1-i}{2}\tau. \quad (A2)$$

We have computed two other solutions, corresponding to $p = 4$ and $5$. These were mentioned in Ref. [20], but not given explicitly. The first one requires $n = 4$ and is given by

$$\begin{aligned}
\tau_1 &= \frac{1}{4} \left( -\frac{1+i}{\sqrt{3}} + 1 - i \right) \tau, \\
\tau_2 &= i \tau_1, \\
\tau_3 &= -i \bar{\tau}_1, \\
\tau_4 &= \bar{\tau}_1.
\end{aligned} \quad (A3)$$

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

$$\begin{aligned}
\tau_1/\tau &= 0.2588533986109182 + 0.0447561340111419i, \\
\tau_2/\tau &= -0.0315468581488038 + 0.2491190542755632i, \\
\tau_3/\tau &= 0.1908290521106672 - 0.2318537492321061i, \\
\tau_4/\tau &= 0.1637288148544367, \\
\tau_5 &= \bar{\tau}_3, \\
\tau_6 &= \bar{\tau}_2, \\
\tau_7 &= \bar{\tau}_1.
\end{aligned} \quad (A4)$$

It should be noted that these solutions are not unique, but we have selected the ones where the error terms, of the order $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 six 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 $O(\tau^p)$ to the error.

<figure><img src="image_17.png" /><figcaption>FIG. 13. Absolute value of the error for the expectation value of the current flowing through the dot, for $N = 6$, $U = 0$, $V = 1.0$, $t = 100$, $J = 0.5$, $J' = 0.5$. The error is computed by comparing the MPS simulations with an exact free-fermion calculation. For small $\tau$ the main source of error is no longer the finite $\tau$, but the successive truncations in the Schmidt decompositions (SVD). This can be seen by comparing the data with cutoff parameter $\delta = 10^{-10}$ to those without any truncation ($\delta = 0$). When the error becomes extremely small, of the order of $10^{-10}$, it stops decreasing with $\tau$. This is due to the finite floating point precision of the machine.</figcaption></figure>

## 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 cutoff 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'$, 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 $\approx 40$.

# APPENDIX B: STEADY STATE IN THE FREE-FERMION CASE

## 1. Steady current

In the free-fermion case ($U = 0$), the transmission and reflexion coefficients $T(k)$ and $R(k)$ for an incident fermion with momentum $k$ have been computed by Branschadel et al. [32]:

$$T(\epsilon) = \frac{1 - [\epsilon/(2J)]^2}{1 + \epsilon^2(J^2 - 2J^2)/(4J^4)}, \quad (B1)$$

$$R(\epsilon) = 1 - T(\epsilon), \quad (B2)$$

$$\epsilon(k) = -2J \cos(k). \quad (B3)$$

Combined with a Landauer approach [33–35] this gives the steady current

$$I = \int_{k_F^-}^{k_F^+} \frac{dk}{2\pi} T(\epsilon(k))v(k), \quad (B4)$$
# APPENDIX A: DETAILS ABOUT THE NUMERICAL SIMULATIONS

## 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 $Δ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 $δ$ 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 $τ = 0.2$ and we used the MPO-based approximation to $\exp(-iHτ)$ that is noted $W^{II}$ in Ref. [20]. As explained in this previous work, one can combine several $W^{II}$ to get a smaller error:

$$W^{II}(\tau_1)W^{II}(\tau_2)\dots W^{II}(\tau_n) = \exp(H\tau) + O(\tau^p). \quad (A1)$$

The simplest case is to use $n = 2$ steps to reduce the error to $O(τ^3)$, and one solution is

$$\tau_1 = \frac{1+i}{2}\tau, \quad \tau_2 = \frac{1-i}{2}\tau. \quad (A2)$$

We have computed two other solutions, corresponding to $p = 4$ and $5$. These were mentioned in Ref. [20], but not given explicitly. The first one requires $n = 4$ and is given by

$$\begin{aligned}
\tau_1 &= \frac{1}{4} \left( -\frac{1+i}{\sqrt{3}} + 1 - i \right) \tau, \\
\tau_2 &= i \tau_1, \\
\tau_3 &= -i \bar{\tau}_1, \\
\tau_4 &= \bar{\tau}_1.
\end{aligned} \quad (A3)$$

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

$$\begin{aligned}
\tau_1/\tau &= 0.2588533986109182 + 0.0447561340111419i, \\
\tau_2/\tau &= -0.0315468581488038 + 0.2491190542755632i, \\
\tau_3/\tau &= 0.1908290521106672 - 0.2318537492321061i, \\
\tau_4/\tau &= 0.1637288148544367, \\
\tau_5 &= \bar{\tau}_3, \\
\tau_6 &= \bar{\tau}_2, \\
\tau_7 &= \bar{\tau}_1.
\end{aligned} \quad (A4)$$

It should be noted that these solutions are not unique, but we have selected the ones where the error terms, of the order $O(τ^p)$, have the smallest prefactors. We have checked the precision of these three different Trotter schemes on a small free-fermion system with six 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 $τ^{p-1}$, due to the fact that the total number of steps is $t/τ$ and each step contributes an amount of order $O(τ^p)$ to the error.

<figure><img src="image_16.png" /><figcaption>FIG. 13. Absolute value of the error for the expectation value of the current flowing through the dot, for $N = 6$, $U = 0$, $V = 1.0$, $t = 100$, $J = 0.5$, $J' = 0.5$. The error is computed by comparing the MPS simulations with an exact free-fermion calculation. For small $τ$ the main source of error is no longer the finite $τ$, but the successive truncations in the Schmidt decompositions (SVD). This can be seen by comparing the data with cutoff parameter $δ = 10^{-10}$ to those without any truncation ($δ = 0$). When the error becomes extremely small, of the order of $10^{-10}$, it stops decreasing with $τ$. This is due to the finite floating point precision of the machine.</figcaption></figure>

## 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 cutoff equal to $δ$ on the discarded weight (the sum of the discarded Schmidt values should be equal to or smaller than $δ$). The effect of the truncation parameter $δ$ is illustrated in Figs. 14 and 15. In this case, a value $δ = 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 $δ$. For smaller $J'$, as shown in Fig. 15, one has to use $δ = 10^{-8}$ in order to get some accurate estimate of the entropy rate. In that case, the bond dimension can exceed 1000 at time $\approx 40$.

# APPENDIX B: STEADY STATE IN THE FREE-FERMION CASE

## 1. Steady current

In the free-fermion case ($U = 0$), the transmission and reflexion coefficients $T(k)$ and $R(k)$ for an incident fermion with momentum $k$ have been computed by Branschadel et al. [32]:

$$T(\epsilon) = \frac{1 - [\epsilon/(2J)]^2}{1 + \epsilon^2(J^2 - 2J'^2)/(4J'^4)}, \quad (B1)$$

$$R(\epsilon) = 1 - T(\epsilon), \quad (B2)$$

$$\epsilon(k) = -2J \cos(k). \quad (B3)$$

Combined with a Landauer approach [33–35] this gives the steady current

$$I = \int_{k_F}^{k_F^+} \frac{dk}{2\pi} T(\epsilon(k))v(k), \quad (B4)$$
<figure><img src="image_18.png" /><figcaption>FIG. 14. Top panel: entanglement entropy of the left lead as a function of time for different values of the SVD truncation parameter $\delta$ which defines the maximum discarded weight at each time step. Middle: maximum MPS bond dimension. Bottom: current flowing through the dot. Parameters of the model: $N = 200$, $U = 2$, $J' = 0.2$, and $V = 1$. The (red) horizontal dotted line indicates the exact value of the steady current at the self-dual point [7]. For these parameters, the upper panel shows that a truncation value as low as $\delta = 10^{-7}$ is required to obtain some accurate value for the entanglement entropy, leading to large bond dimensions (above 1000). In contrast, $\delta = 10^{-6}$ appears to be enough to estimate the current correctly. Simulations performed with Trotter step $\tau = 0.2$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 15. Same as Fig. 14, for $J' = 0.05$. For these parameters, a truncation value as low as $\delta = 10^{-8}$ is required to obtain some accurate value for the entanglement entropy, up to $t = 60$. One observes some small amplitude oscillations in $S(t)$, with the same period $T_{osc} = 4\pi/V$, as for $I(t)$.</figcaption></figure>

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). \qquad (B5)$$

For the IRLM, using Eq. (B1), the integral gives

$$2\pi I = -V \frac{J'^4}{a} + 4J'^2 \frac{(1 - J'^2)^2}{a^{\frac{3}{2}}} \arctan \left( \frac{\sqrt{a}}{4J'^2} V \right), \\ a = 1 - 2J'^2, \qquad (B6)$$

where we assumed $J = 1$, $V < 4$, and $J'^2 < 1$. This function is plotted in Fig. 16 for three different values of $J'$. In

the (scaling) limit where $J' \ll 1$ we define $T_B = (J')^2$ and Eq. (B7) simplifies to

$$\frac{2\pi I}{T_B} = 4 \arctan \left( \frac{V}{4T_B} \right), \qquad (B7)$$

in agreement with Ref. [7]. Note that the expression $I(t)$ of the current at finite time and $U = 0$ can be found in the Appendix of Ref. [30]. It shows damped oscillations at frequency $T_{osc} = 4\pi/V$, a relaxation time $O(J'^2)$, and converges to Eq. (B7) when $t \to \infty$.

## 2. Density drop across the dot

We make the approximation that the fermions are pointlike particles which propagate ballistically in leads, at some group velocity $v(k)$. This is a semiclassical approximation (called hydrodynamical approximation in Ref. [36]) 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
<figure><img src="image_17.png" /><figcaption>FIG. 14. Top panel: entanglement entropy of the left lead as a function of time for different values of the SVD truncation parameter $\delta$ which defines the maximum discarded weight at each time step. Middle: maximum MPS bond dimension. Bottom: current flowing through the dot. Parameters of the model: $N = 200$, $U = 2$, $J' = 0.2$, and $V = 1$. The (red) horizontal dotted line indicates the exact value of the steady current at the self-dual point [7]. For these parameters, the upper panel shows that a truncation value as low as $\delta = 10^{-7}$ is required to obtain some accurate value for the entanglement entropy, leading to large bond dimensions (above 1000). In contrast, $\delta = 10^{-6}$ appears to be enough to estimate the current correctly. Simulations performed with Trotter step $\tau = 0.2$.</figcaption></figure>

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). \qquad (B5)$$

For the IRLM, using Eq. (B1), the integral gives

$$2\pi I = -V \frac{J'^4}{a} + 4J'^2 \frac{(1 - J'^2)^2}{a^{\frac{3}{2}}} \arctan \left( \frac{\sqrt{a}}{4J'^2} V \right), \\ a = 1 - 2J'^2, \qquad (B6)$$

where we assumed $J = 1$, $V < 4$, and $J'^2 < 1$. This function is plotted in Fig. 16 for three different values of $J'$. In

<figure><img src="image_18.png" /><figcaption>FIG. 15. Same as Fig. 14, for $J' = 0.05$. For these parameters, a truncation value as low as $\delta = 10^{-8}$ is required to obtain some accurate value for the entanglement entropy, up to $t = 60$. One observes some small amplitude oscillations in $S(t)$, with the same period $T_{osc} = 4\pi/V$, as for $I(t)$.</figcaption></figure>

the (scaling) limit where $J' \ll 1$ we define $T_B = (J')^2$ and Eq. (B7) simplifies to

$$\frac{2\pi I}{T_B} = 4 \arctan \left( \frac{V}{4T_B} \right), \qquad (B7)$$

in agreement with Ref. [7]. Note that the expression $I(t)$ of the current at finite time and $U = 0$ can be found in the Appendix of Ref. [30]. It shows damped oscillations at frequency $T_{osc} = 4\pi/V$, a relaxation time $O(J'^2)$, and converges to Eq. (B7) when $t \to \infty$.

## 2. Density drop across the dot

We make the approximation that the fermions are pointlike particles which propagate ballistically in leads, at some group velocity $v(k)$. This is a semiclassical approximation (called hydrodynamical approximation in Ref. [36]) 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
<figure><img src="image_20.png" /><figcaption>FIG. 16. Rescaled steady current as a function of the rescaled voltage in the free case ($U=0$), using $T_B = J'^2$. The symbols represent the DMRG data for three different values of $J'$, and the full lines are, respectively, the exact result [Eq. (B7)] for $J' = 0.5$ (green), for $J' = 0.35$ (blue), for $J' = 0.2$ (orange), and the limit when $J'$ is small [Eq. (B7), black]. The DMRG data are in perfect agreement with the exact results, but since chosen $J'$ are not very small one observes some deviations from Eq. (B7).</figcaption></figure>

momentum distributions and the scattering on the dot, we obtain

$$n(k)^L = \begin{cases} 0, & \text{if } k_F^+ < k \\ 1, & \text{if } -k_F^- < k < k_F^+ \\ \mathcal{R}(k), & \text{if } -k_F^+ < k < -k_F^- \\ 0, & \text{if } k < -k_F^+ \end{cases} \quad (B8)$$

$$n(k)^R = \begin{cases} 0, & \text{if } k_F^+ < k \\ \mathcal{T}(k), & \text{if } -k_F^- < k < k_F^+ \\ 1, & \text{if } -k_F^- < k < -k_F^- \\ 0, & \text{if } k < -k_F^- \end{cases} \quad (B9)$$

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

$$\rho^L = \frac{1}{2} + \int_{k_F^-}^{k_F^+} \frac{dk}{2\pi} \mathcal{R}(\epsilon(k)), \quad (B10)$$

$$\rho^R = \frac{1}{2} - \int_{k_F^-}^{k_F^+} \frac{dk}{2\pi} \mathcal{R}(\epsilon(k)). \quad (B11)$$

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 quasisteadyregion 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)) \quad (B12)$$

across the dot.

### 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 [\mathcal{T}(\epsilon) \ln \mathcal{T}(\epsilon) + \mathcal{R}(\epsilon) \ln \mathcal{R}(\epsilon)]. \quad (B13)$$

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

In the scaling regime $J'$ is small, the energy $\epsilon$ are small (because the bias $V$ is small), but the ratio $\epsilon/J'^2$ is of order one. In this limit, the transmission coefficient [Eq. (B1)] becomes

$$\mathcal{T}(x) = \frac{1}{1 + x^2}, \quad (B14)$$

where $x = \epsilon/(2T_B)$ and $T_B = J'^2$. With a change of variable, the entropy rate [Eq. (B13)] can be expressed as

$$\alpha(v) = \frac{2T_B}{2\pi} \int_{-v/4}^{v/4} dx \left[ \frac{1}{1+x^2} \ln(1+x^2) + \frac{x^2}{1+x^2} \ln \left( \frac{1+x^2}{x^2} \right) \right], \quad (B15)$$

where $v = V/T_B$. The integral can be computed explicitly and the final result is

$$\frac{\alpha(v)}{T_B} = \frac{2}{\pi} \left[ 2[1 + \ln(v/4)] \arctan(v/4) + \frac{1}{4} v \ln(v^2 + 16) - \frac{1}{2} v \ln(v) - i \mathrm{Li}_2\left(-\frac{i v}{4}\right) + i \mathrm{Li}_2\left(\frac{i v}{4}\right) \right], \quad (B16)$$

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

$$\alpha(v \to \infty)/T_B = 2. \quad (B17)$$

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 + O(v^5). \quad (B18)$$

[1] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, *Rev. Mod. Phys.* **83**, 863 (2011).

[2] J. Eisert, M. Friesdorf, and C. Gogolin, *Nat. Phys.* **11**, 124 (2011).

[3] U. Schollwöck, *Ann. Phys.* **326**, 96 (2011).

[4] A. C. Hewson, *The Kondo Problem to Heavy Fermions* (Cambridge University Press, Cambridge, UK, 1993).
<figure><img src="image_19.png" /><figcaption>FIG. 16. Rescaled steady current as a function of the rescaled voltage in the free case ($U=0$), using $T_B = J'^2$. The symbols represent the DMRG data for three different values of $J'$, and the full lines are, respectively, the exact result [Eq. (B7)] for $J' = 0.5$ (green), for $J' = 0.35$ (blue), for $J' = 0.2$ (orange), and the limit when $J'$ is small [Eq. (B7), black]. The DMRG data are in perfect agreement with the exact results, but since chosen $J'$ are not very small one observes some deviations from Eq. (B7).</figcaption></figure>

momentum distributions and the scattering on the dot, we obtain

$$n(k)^L = \begin{cases} 0, & \text{if } k_F^+ < k \\ 1, & \text{if } -k_F^- < k < k_F^+ \\ \mathcal{R}(k), & \text{if } -k_F^+ < k < -k_F^- \\ 0, & \text{if } k < -k_F^+ \end{cases} \quad (B8)$$

$$n(k)^R = \begin{cases} 0, & \text{if } k_F^+ < k \\ \mathcal{T}(k), & \text{if } -k_F^- < k < k_F^+ \\ 1, & \text{if } -k_F^- < k < -k_F^- \\ 0, & \text{if } k < -k_F^- \end{cases} \quad (B9)$$

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

$$\rho^L = \frac{1}{2} + \int_{k_F^-}^{k_F^+} \frac{dk}{2\pi} \mathcal{R}(\epsilon(k)), \quad (B10)$$

$$\rho^R = \frac{1}{2} - \int_{k_F^-}^{k_F^+} \frac{dk}{2\pi} \mathcal{R}(\epsilon(k)). \quad (B11)$$

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 quasisteadyregion 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)) \quad (B12)$$

across the dot.

### 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 [\mathcal{T}(\epsilon) \ln \mathcal{T}(\epsilon) + \mathcal{R}(\epsilon) \ln \mathcal{R}(\epsilon)]. \quad (B13)$$

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

In the scaling regime $J'$ is small, the energy $\epsilon$ are small (because the bias $V$ is small), but the ratio $\epsilon/J'^2$ is of order one. In this limit, the transmission coefficient [Eq. (B1)] becomes

$$\mathcal{T}(x) = \frac{1}{1 + x^2}, \quad (B14)$$

where $x = \epsilon/(2T_B)$ and $T_B = J'^2$. With a change of variable, the entropy rate [Eq. (B13)] can be expressed as

$$\alpha(v) = \frac{2T_B}{2\pi} \int_{-v/4}^{v/4} dx \left[ \frac{1}{1+x^2} \ln(1+x^2) + \frac{x^2}{1+x^2} \ln \left( \frac{1+x^2}{x^2} \right) \right], \quad (B15)$$

where $v = V/T_B$. The integral can be computed explicitly and the final result is

$$\frac{\alpha(v)}{T_B} = \frac{2}{\pi} \left[ 2[1 + \ln(v/4)] \arctan(v/4) + \frac{1}{4} v \ln(v^2 + 16) - \frac{1}{2} v \ln(v) - i \mathrm{Li}_2\left(-\frac{iv}{4}\right) + i \mathrm{Li}_2\left(\frac{iv}{4}\right) \right], \quad (B16)$$

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

$$\alpha(v \to \infty)/T_B = 2. \quad (B17)$$

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 + O(v^5). \quad (B18)$$

[1] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, *Rev. Mod. Phys.* **83**, 863 (2011).

[2] J. Eisert, M. Friesdorf, and C. Gogolin, *Nat. Phys.* **11**, 124 (2011).

[3] U. Schollwöck, *Ann. Phys.* **326**, 96 (2011).

[4] A. C. Hewson, *The Kondo Problem to Heavy Fermions* (Cambridge University Press, Cambridge, UK, 1993).
[5] P. B. Wiegmann and A. M. Finkel'shtein, Zh. Eksp. Teor. Fiz. **75**, 204 (1978) [Sov. Phys.–JETP **48**, 102 (1978)].

[6] V. M. Filyov and P. B. Wiegmann, Phys. Lett. A **76**, 283 (1980).

[7] E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. Lett. **101**, 140601 (2008).

[8] P. Fendley, A. W. W. Ludwig, and H. Saleur, Phys. Rev. Lett. **74**, 3005 (1995).

[9] P. Fendley, A. W. W. Ludwig, and H. Saleur, Phys. Rev. B **52**, 8934 (1995).

[10] P. Mehta and N. Andrei, Phys. Rev. Lett. **96**, 216802 (2006).

[11] D. Gobert, C. Kollath, U. Schollwöck, and G. Schütz, Phys. Rev. E **71**, 036102 (2005).

[12] T. Sabetta and G. Misguich, Phys. Rev. B **88**, 245114 (2013).

[13] P. Schmitteckert, Phys. Rev. B **70**, 121302 (2004).

[14] G. Schneider and P. Schmitteckert, arXiv:cond-mat/0601389.

[15] G. Vidal, Phys. Rev. Lett. **93**, 040502 (2004).

[16] S. R. White and A. E. Feiguin, Phys. Rev. Lett. **93**, 076401 (2004).

[17] A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, J. Stat. Mech. (2004) P04005.

[18] N. Laflorencie, Phys. Rep. **646**, 1 (2016).

[19] ITensor Library, http://itensor.org (version 2.0)

[20] M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore, and F. Pollmann, Phys. Rev. B **91**, 165112 (2015).

[21] V. Eisler and I. Peschel, Europhys. Lett. **99**, 20001 (2012).

[22] D. M. Kennes, V. Meden, and R. Vasseur, Phys. Rev. B **90**, 115101 (2014).

[23] R. Vasseur and H. Saleur, Sci. Post Phys. **3**, 001 (2017).

[24] P. Calabrese and J. Cardy, J. Stat. Mech. (2007) P10004.

[25] J.-M. Stéphan and J. Dubail, J. Stat. Mech. (2011) P08019.

[26] I. Klich and L. Levitov, Phys. Rev. Lett. **102**, 100502 (2009).

[27] H. F. Song, S. Rachel, C. Flindt, I. Klich, N. Laflorencie, and K. Le Hur, Phys. Rev. B **85**, 035409 (2012).

[28] H. F. Song, C. Flindt, S. Rachel, I. Klich, and K. Le Hur, Phys. Rev. B **83**, 161408 (2011).

[29] C. Karrasch, M. Pletyukhov, L. Borda, and V. Meden, Phys. Rev. B **81**, 125122 (2010).

[30] Y. Vinkler-Aviv, A. Schiller, and F. B. Anders, Phys. Rev. B **90**, 155110 (2014).

[31] L. Freton, E. Boulat, and H. Saleur, Nucl. Phys. B **874**, 279 (2013).

[32] A. Branschädel, E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. B **82**, 205414 (2010).

[33] R. Landauer, IBM J. Res. Dev. **1**, 223 (1957).

[34] Y. Imry and R. Landauer, Rev. Mod. Phys. **71**, S306 (1999).

[35] Y. V. Nazarov and Y. Blanter, *Quantum Transport: Introduction to Nanoscience* (Cambridge University Press, Cambridge, UK, 2009).

[36] T. Antal, P. L. Krapivsky, and A. Rákos, Phys. Rev. E **78**, 061115 (2008).
[5] P. B. Wiegmann and A. M. Finkel'shtein, Zh. Eksp. Teor. Fiz. **75**, 204 (1978) [Sov. Phys.–JETP **48**, 102 (1978)].

[6] V. M. Filyov and P. B. Wiegmann, Phys. Lett. A **76**, 283 (1980).

[7] E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. Lett. **101**, 140601 (2008).

[8] P. Fendley, A. W. W. Ludwig, and H. Saleur, Phys. Rev. Lett. **74**, 3005 (1995).

[9] P. Fendley, A. W. W. Ludwig, and H. Saleur, Phys. Rev. B **52**, 8934 (1995).

[10] P. Mehta and N. Andrei, Phys. Rev. Lett. **96**, 216802 (2006).

[11] D. Gobert, C. Kollath, U. Schollwöck, and G. Schütz, Phys. Rev. E **71**, 036102 (2005).

[12] T. Sabetta and G. Misguich, Phys. Rev. B **88**, 245114 (2013).

[13] P. Schmitteckert, Phys. Rev. B **70**, 121302 (2004).

[14] G. Schneider and P. Schmitteckert, arXiv:cond-mat/0601389.

[15] G. Vidal, Phys. Rev. Lett. **93**, 040502 (2004).

[16] S. R. White and A. E. Feiguin, Phys. Rev. Lett. **93**, 076401 (2004).

[17] A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, J. Stat. Mech. (2004) P04005.

[18] N. Laflorencie, Phys. Rep. **646**, 1 (2016).

[19] ITensor Library, http://itensor.org (version 2.0)

[20] M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore, and F. Pollmann, Phys. Rev. B **91**, 165112 (2015).

[21] V. Eisler and I. Peschel, Europhys. Lett. **99**, 20001 (2012).

[22] D. M. Kennes, V. Meden, and R. Vasseur, Phys. Rev. B **90**, 115101 (2014).

[23] R. Vasseur and H. Saleur, Sci. Post Phys. **3**, 001 (2017).

[24] P. Calabrese and J. Cardy, J. Stat. Mech. (2007) P10004.

[25] J.-M. Stéphan and J. Dubail, J. Stat. Mech. (2011) P08019.

[26] I. Klich and L. Levitov, Phys. Rev. Lett. **102**, 100502 (2009).

[27] H. F. Song, S. Rachel, C. Flindt, I. Klich, N. Laflorencie, and K. Le Hur, Phys. Rev. B **85**, 035409 (2012).

[28] H. F. Song, C. Flindt, S. Rachel, I. Klich, and K. Le Hur, Phys. Rev. B **83**, 161408 (2011).

[29] C. Karrasch, M. Pletyukhov, L. Borda, and V. Meden, Phys. Rev. B **81**, 125122 (2010).

[30] Y. Vinkler-Aviv, A. Schiller, and F. B. Anders, Phys. Rev. B **90**, 155110 (2014).

[31] L. Freton, E. Boulat, and H. Saleur, Nucl. Phys. B **874**, 279 (2013).

[32] A. Branschädel, E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. B **82**, 205414 (2010).

[33] R. Landauer, IBM J. Res. Dev. **1**, 223 (1957).

[34] Y. Imry and R. Landauer, Rev. Mod. Phys. **71**, S306 (1999).

[35] Y. V. Nazarov and Y. Blanter, *Quantum Transport: Introduction to Nanoscience* (Cambridge University Press, Cambridge, UK, 2009).

[36] T. Antal, P. L. Krapivsky, and A. Rákos, Phys. Rev. E **78**, 061115 (2008).