{'fraction_non_alphanumeric': 0.03163652081393802, 'fraction_numerical': 0.010273972602739725, 'mean_word_length': 4.470443797744634, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical method for non-linear steady-state transport in one-dimensional correlated conductors

We present a method for investigating the steady-state transport properties of one-dimensional correlated quantum systems. Using a procedure based on our analysis of finite-size effects in a related classical model (LC line) we show that stationary currents can be obtained from transient currents in finite systems driven out of equilibrium. The non-equilibrium dynamics of correlated quantum systems is calculated using the time-evolving block decimation method. To demonstrate our method we determine the full I-V characteristic of the spinless fermion model with nearest-neighbour hopping t_H and interaction V_H using two different setups to generate currents (turning on/off a potential bias). Our numerical results agree with exact results for non-interacting fermions (V_H=0). For interacting fermions we find that in the linear regime eV<<4t_H the current I is independent from the setup and our numerical data agree with the predictions of the Luttinger liquid theory combined with the Bethe Ansatz solution. For larger potentials V the steady-state current depends on the current-generating setup and as V increases we find a negative differential conductance with one setup while the currents saturate at finite values in the other one. Both effects are due to finite renormalized bandwidths.

Numerical method for non-linear steady-state transport in one-dimensional correlated conductors We present a method for investigating the steady-state transport properties of one-dimensional correlated quantum systems. Using a procedure based on our analysis of finite-size effects in a related classical model (LC line) we show that stationary currents can be obtained from transient currents in finite systems driven out of equilibrium. The non-equilibrium dynamics of correlated quantum systems is calculated using the time-evolving block decimation method. To demonstrate our method we determine the full I-V characteristic of the spinless fermion model with nearestneighbour hopping tH and interaction VH using two different setups to generate currents (turning on/off a potential bias). Our numerical results agree with exact results for non-interacting fermions (VH = 0). For interacting fermions we find that in the linear regime eV ≪ 4tH the current I is independent from the setup and our numerical data agree with the predictions of the Luttinger liquid theory combined with the Bethe Ansatz solution. For larger potentials V the steady-state current depends on the current-generating setup and as V increases we find a negative differential conductance with one setup while the currents saturate at finite values in the other one. Both effects are due to finite renormalized bandwidths.

I. INTRODUCTION
Much of our understanding of electronic transport in solids is based on a picture of weakly-interacting charge carriers such as Fermi liquid quasi-particles. Onedimensional electron systems are well-known examples where this approximation fails. The low-energy physics of these systems is described by the Tomonaga-Luttinger liquid (TLL) theory and is realized in quantum wires such as carbon nanotubes, nanowires in semiconductor heterostructures, or metal atomic chains on surface substrates. [1][2][3] The transport properties of one-dimensional correlated systems have been extensively studied during the last two decades. [2][3][4][5][6] A major goal is to determine and understand the current-voltage characteristic I − V of various systems made of quantum dots and wires. Most of these studies have been restricted to a regime where the energy scale of the current-generating electromagnetic field or potential bias is small in comparison to the energy scale of the unperturbated systems (i.e., the band width in lattice models and the Fermi velocity in continuum models). The non-linear regime has mostly been investigated in quantum contact problems, where the interaction is confined to a small region of the system. For instance, non-linear current-voltage characteristics have been calculated within the TLL theory such as the powerlaw I ∼ V α behaviour for the transport through a weak link 7,8 or the exact I(V ) curve for the current through a point contact in a fractional quantum Hall edge state device. 9,10 However, the TLL theory is limited to lowenergy excitations with linear dispersion and thus to potential biases V which are weak compared to the band width. Only recently the implications of a non-linear dispersion have started to be considered. 11 Thus current-voltage characteristics from the TLL theory are actually limited to a weak-bias regime.
There are few works presenting full current-voltage characteristics with a voltage up to the largest energy scale of the system (for instance, see Refs. 12, 13, and 14) and none is concerned with one-dimensional correlated conductors. Thus transport properties beyond linear response are poorly understood. We believe that it is important to attain a better knowledge of the non-linear transport properties in one-dimensional correlated conductors. First, the validity of weak-bias approaches can be confirmed only if one obtains some quantitative estimates of non-linear effects. Moreover, non-linear devices play a significant role in electronics and studies of nonlinear dynamics are required to reveal the full potential functionality of quantum wires as electronic circuit components.
In this work we develop and apply a method for investigating the zero-temperature DC transport properties of one-dimensional correlated conductors for potential biases up to the order of the band width. For this purpose we study a well-known one-dimensional lattice model described by the half-filled spinless fermion Hamiltonian with nearest-neighbour repulsion. Even though this model is exactly solvable by the Bethe Ansatz and the low-energy physics are determined by the generic TLL phenomenology 2,3 , its transport properties in the non-linear regime cannot be obtained analytically.
To determine the transport properties we simulate the quantum dynamics of single chains which are driven out of equilibrium by a potential bias between the left-and right-hand halves of the chain and calculate the resulting currents through the middle bond. The simulation of out-of-equilibrium quantum many-particle systems is one of the major challenges in computational physics.
Recently, a family of numerical methods has been developed to simulate the real-time evolution of quantum lattice systems such as one-dimensional Hamiltonians with short-range interactions. [15][16][17] The most prominent ones are the time-dependent Density-Matrix Renormalization Group (td-DMRG) and the time-evolving block decimation (TEBD) method. 18,19 Various flavours of td-DMRG have been successfully applied in studies of quantum many-body dynamics. In particular, they have proven to be promising tools for investigating electronic transport in strongly correlated nanostructures and onedimensional conductors. 12,13,[20][21][22][23][24][25][26][27][28][29][30] Surprisingly, the original TEBD method 18,19 has not been applied to electronic transport problems yet. Both td-DMRG and TEBD methods can be described within a common mathematical framework, the matrix product quantum states. 17,31,32 Their accuracy and efficiency depend essentially on the amount of entanglement in the quantum system and should be similar. Admittedly, the TEBD algorithm is restricted to a small family of systems (one-dimensional Hamiltonians or ladder systems with nearest-neighbour interactions only) while td-DMRG techniques are more versatile. However, the TEBD algorithm is naturally parallelizable and thus fully scalable, which is a required feature in high-performance computing, while the efficient parallelization of the DMRG algorithms remains an open challenge. 33,34 In this work we employ the TEBD method for computing the non-equilibrium dynamics of the spinless fermion model.
With the TEBD method we can compute the non equilibrium quantum dynamics of lattice models with a finite number of sites N over a finite period of time t. As we are primarily interested in determining the DC transport properties, TEBD results must be extrapolated to the thermodynamic limit N → ∞ and to the steady-state limit t → ∞. However, finite-size and finite-time effects are very complex in these out-of-equilibrium quantum systems (for instance, see Ref. 35) and extrapolations are difficult because we know very little about the scaling of currents with system size N and time t. For this reason we have investigated this scaling in an exactly solvable classical model, the so-called LC line. Using this information we have developed a method for quantum systems which allows us to obtain reliable quantitative results for stationary currents I from numerical data for rather small system sizes and short simulation times.
In this paper we show that our extrapolation approach allows one to determine the full I−V curves of interacting one-dimensional conductors using the spinless fermion model for illustration. In the non-interacting case this model can be solved exactly using the equation of motion method. The outcomes of this special case confirm our extrapolation method and reveal a negligible numerical error for the TEBD simulation results. For interacting fermions comparisons with predictions from the TLL theory and the Bethe-Ansatz solution confirm the validity of our method in the linear response regime. While for the linear regime the specific setup does not matter, it is highly decisive for the non-linear current-voltage characteristics. We basically distinguish between two different ways of creating a current flow. 35 In the first setup (I) we apply an initial voltage and calculate the ground state which has different particle numbers in its two halves and then let the system evolve with an overall equal onsite potential. In the second setup (II) we calculate the ground state without a potential difference but turn it on for the real time evolution. We have found that while the general shape of the current-curve as a function of time is dominated by classical effects (i.e., which are also found in the LC line model), the current-voltage characteristics are primarily determined by the chosen setup. For setup (I) the system shows a positive differential conductance for the full voltage-range and saturates at a finite value for very large potential differences. For the second setup (II) the linear response coincides with the one for setup (I) but for higher voltages we observe a negative differential conductance. Both effects are also present in the non-interacting case and come from the finite bandwidth and the non-linear dispersion (i.e., the energy-dependent density of states) of the excitations. 35,36 Our paper is organized as follows: In the next section we introduce the spinless fermion model, and the methods used in this work. In the third section we investigate the classical LC line to explain the basic behaviour of the current in finite and infinite one-dimensional chains and derive our method to extrapolate the stationary current from finite system results. In the fourth section we describe finite-size effects and the convergence to a stationary current for infinite system sizes, while the I − V characteristic of the spinless fermion model and comparisons with exact results are shown in the fifth section. Finally, we summarize our findings in the last section. Some calculations are detailed in appendices.

A. Spinless fermion model
We consider a one-dimensional lattice model representing correlated conductors driven out of equilibrium by a potential bias. For spinless fermions the Hamiltonian without a potential bias is where t H denotes the hopping amplitude between nearest-neighbour sites, V H is the Coulomb repulsion between spinless fermions on nearest-neighbour sites, and n j = c † j c j . At half filling (N/2 fermions in the N -site lattice) this Hamiltonian describes an ideal conductor for One-dimensional conductor consisting of two coupled leads. A bias between its right-and left-hand halves is applied and the current is measured in the middle of the system at the junction between both sides (dashed connection).
−2t H < V H ≤ 2t H . It can be interpreted as a system of spin-polarized electrons. The low-energy behaviour of this lattice model is described by the TLL theory. 2,3 The generic properties of a TLL are determined by two parameters: The velocity of elementary excitations (renormalized Fermi velocity) v and a dimensionless parameter K. From the Bethe Ansatz solution we know the relation between TLL parameters and the parameters t H and V H of the spinless fermion model at half filling and where a L is the lattice constant. To drive the system out of equilibrium we a use step-like potential bias between the left-and right-hand halves of the chain The potential energy step is set by ∆ǫ = |eV | where V is the voltage bias and e is the elementary charge. It is possible to use a smoother potential profile but the results for the stationary current using our extrapolation method are only slightly affected by the specific shape as long as the non-constant part in the middle is rather smooth and locally confined.
The above system can be seen as two coupled interacting leads made of the sites {1 . . . N 2 } and { N 2 + 1 . . . N }, respectively, see Fig. 1. The coupling is given by a hopping term t ′ H between the left-and right-hand sides of the system (i.e., between site N/2 and site N/2 + 1) and an additional coupling V ′ H ( n N/2 − 1/2)( n N/2+1 − 1/2). Here we discuss the homogeneous system only (t ′ H = t H and V ′ H = V H ) but our approach can be easily extended to systems with a weak link t ′ H < t H representing the tunneling through a nanostructure 20 or with a site representing a quantum dot, such as the Interacting Resonant Level Model (IRLM) 12 , as well as to systems with a few additional sites intercalated between both leads and representing nanostructures with internal degrees of freedom. 23 The current operator between the pair of sites (k, k+1) is For a given (time-dependent) quantum state we define the current flowing between both halves of the system (see Fig. 1) as the expectation value of the current operator for the site pair in the middle of the system We note that where n k (8) are the (time-dependent) charges in the left-and righthand halves of the chain, respectively. As the number of particles n is conserved in our models, Q L (t) + Q R (t) = −e · n = const. The stationary current is a constantly flowing current in an infinitely large system after the settling timeJ = lim t→∞ lim N →∞ J(t) (9) and according to the TLL theory in the linear regime (small ∆ǫ) it is given by 2 We set t H = a L = 1 for all numerical simulations and, if units are not given explicitly, e =h = 1.
In this work both setups used to generate a current (see the next section) lead to an overall half-filled system. Thus in the weak potential bias regime both system halves remain approximately half-filled at all times. Consequently, we expect that the current is given by the Luttinger liquid prediction (10) together with the Bethe Ansatz result (3) in the linear response regime.

B. Setups for non-equilibrium simulations
We employ two different setups to generate currents in the lattice models. 35 In the first one (I) we The charge reservoirs (two halves: source and drain) have different potentials but are coupled for t = 0 while for t > 0 the potential difference is set to ∆ǫ = 0 instantaneously.
according to the Hamiltonian H 0 (i.e. without potential bias) This setup describes an inhomogeneous initial state with more particles in one half of the system than in the other one. Thus particles flow from one side to the other one for t > 0. It corresponds to a one-dimensional scattering experiment in which particles are emitted on one side of the system with energies between [−∆ǫ/2, ∆ǫ/2], scattered at the junction between both system halves, and then (partially) transmitted to the opposite side. This picture of transport through junctions is often used in theoretical investigations.
In the second setup (II) we prepare the system at time t = 0 in the ground state |φ(∆ǫ = 0) of the Hamiltonian H 0 (i.e., without potential bias). For later times t > 0 the time evolution of the system is determined by the Hamiltonian H = H 0 − H B , i.e. with a potential bias that causes the current to flow in the same direction as in setup (I) see Fig. 3. Thus it corresponds more closely to the actual experimental situation with a voltage source generating a current in a conducting wire. We should point out that in setup (I) the two leads are decoupled with respect to the Coulomb-interaction V ′ H for the calculation of the ground state. Our tests have revealed that otherwise a strong dependency of the stationary current on the system size appears, that is, for smaller system sizes the stationary current becomes higher but for N → ∞ it approaches the same constant value as the value one gets with V ′ H = 0 for t = 0. Choosing V ′ H = 0 for the computation of the ground state therefore decreases the finite-size error.
Generally, currents (6) calculated with the states (11) and (12) are different. In the strong-bias limit |∆ǫ| ≫ t H , V H it is easy to show that the steady-state current remains finite for the first setup while it vanishes for the second one. Recently, it has been reported that initial conditions (quenching an interaction term or a tunneling term) can also alter the steady-state current flowing through a quantum point contact between two TLL leads which have been driven out of equilibrium by an external bias. 37 In the weak-bias limit |∆ǫ| ≪ t H , however, a simple perturbation calculation shows that both setups yield the same linear response for the stationary current. Thus in this regime both setups can be used indifferently but in the non-linear regime we must distinguish them.
In most theoretical studies the potential bias switching is not instantaneous but adiabatic. This can also be used in numerical simulations, for instance, see Refs. 20 and 27. Our tests do not reveal any significant differences for the steady-state current depending on the switching rate as long as the potential changes in a time scale which is much smaller than the time scale associated with the motion of particles from one reservoir to the other one. Since our extrapolation method is designed to work with an instantaneous change in the potential difference and as numerical simulations are simpler with it, we prefer this approach.

C. One-particle equation of motion
Without Coulomb-interaction V H (= V ′ H ) = 0 the spinless fermion model, described by the Hamiltonian (1), reduces to the tight-binding model (non-interacting fermions without spin degree of freedom). A onedimensional chain in the tight-binding model can be described by the single particle reduced density matrix where the time evolution is given by the one-particle equation of motion H (1) t>0 denotes the single particle Hamilton matrix of size N × N with which the system is evolved in time. More precisely H (1) t>0 is the one-particle representation of H 0 or H 0 + H B depending on the specific setup. The particle number expectation values coincide with the diagonal terms of the reduced density matrix (15) and the expectation value for the current operator (5) from site k to k + 1 can be taken from off-diagonal entries Consequently, the dynamics of a tight-binding chain can be computed numerically without additional truncation error with a runtime of O(N 3 ) for any point in time.

D. TEBD method
The TEBD method 18,19 is based on a specific representation of quantum states by matrix product states (MPS). For an N -site lattice it is generally written Sums over an index j k run over a complete basis of the site k (for instance 2 states for the spinless fermion model). Every quantum state of the Fock space associated with a finite lattice can be represented exactly in this form if the matrix dimensions can be as large as the square root of the Fock space dimension (for instance 2 N/2 for the spinless fermion model). In numerical computations, however, the matrix dimension must be kept smaller than a relatively small upper limit χ c . Fortunately, for many one-dimensional systems this truncation is possible and can lead to a dramatic computational speedup while keeping the error in computed observables conveniently low. This is done by using the Schmidt decomposition of the density matrix at each bond to calculate the matrices λ [k] from its eigenvalues and the Γ [k]j k from its eigenvectors. Given a bond k in the chain, the Schmidt decomposition at this bipartite split is defined by The The λ α k are the matrix elements (λ [k] ) α k and either of the |Φ α k can be used to build the matrices Γ [k]j k from equation (17). Thus, it is possible to keep only the largest χ c eigenvalues and throw away the rest while re-normalizing the state, such that the sum of the discarded values is smaller than an arbitrary error ǫ The best candidates are states for which the Schmidt dimension, and hence the dimension of the MPS matrices is low (like product states), or states for which the Schmidt coefficients display an exponential decay, where a large number of eigenvalues can be discarded without significant information loss. The great advantage of the TEBD algorithm is the possibility to compute the time-evolution of a state using a time-dependent Hamiltonian In a numerical implementation, δt has to be discrete, such that the total simulation time τ = n t · δt, where n t is the number of time steps and δt the numerical time step. The time-evolution can be used as well to calculate the ground state |ψ gr of a Hamiltonian H. This is done by taking the time to be imaginary and projecting (effectively 'cooling') a starting state |ψ P to the ground state of H Given a one-dimensional system of size N with nearestneighbour interaction we can split the Hamiltonian into two sums H N = F + G over even and odd sites, with By using the Suzuki-Trotter decomposition for exponential operators we can reduce the process of time-evolution to the successive application of operators which act only on two sites The TEBD algorithm then proceeds to update the MPS representation by calculating the new Schmidt decomposition at the corresponding bond each time after a twosite operator is applied to it. This succession of two-site operators imposes a linear dependence of the computational cost on the system size N . Each time a Schmidt decomposition is computed, we can truncate the dimension of the MPS by keeping only a maximal number of eigenvalues, using a suitable limitation criterion. This is desirable, as the computational cost of both Schmidt decomposition and update of the MPS representation is O(χ c 3 ) as explicitly shown in Refs. 18 and 19. The alternate updating of even and odd sites according to (26) makes the TEBD algorithm highly parallelizable using up to (N − 1)/2 threads with a very low communication overhead.
One obvious error source is the one stemming from the discretization of time in order to numerically compute equation (21), so it is necessary to keep the actual time step δt small enough to have a good approximation of H(t) in the interval [t, t + δt]. But this error source is not inherent to the TEBD implementation. The main error sources in the algorithm are the Suzuki-Trotter approximation and the Schmidt truncation. In order to improve the time step error, one can use higher dimensional Suzuki-Trotter formulas, at a computational cost which scales linearly with higher-order approximations 38 , or we can decrease the time step δt, which also comes at a linear cost n t = δt τ . The dominating error in probably all setups of interest is thus the truncation error. It is also very difficult to compensate for this error, because of the O(χ c 3 ) scaling of the computational cost. As a trivial example, when starting a real time evolution with a state with Schmidt dimension smaller than χ c , the simulation runs with a constant truncation error. As the simulation continues, there comes a point where more states than χ c are needed, and the truncation error quickly overcomes the Trotter error, which basically defines a runaway time for the simulation, as reported in Ref. 39. A reasonable and often used estimator of the truncation error is the discarded weight where D k denotes the discarded weight for a bipartite split at the k-th bond and χ m < χ c is the number of kept eigenvalues. For our real time simulations, we use a maximal Schmidt dimension in the range of 300 ≤ χ c ≤ 500, depending on the specific parameter combination. We use a site-dependent χ c (site) and allow our simulations to adapt χ c (site) if test values (discarded weight, von-Neumann entropy) show the necessity to do so. We have to remark that for higher |V H | and ∆ǫ the correlations (and thus the needed Schmidt dimension) within the spinless fermion model grow very quickly. Hence, the chosen χ c -range can lead to significantly larger errors in the regime |V H | ≥ 1.6t H and ∆ǫ ≥ 3t H . Figure 4 shows the discarded weight for different parameters for a split in the middle of the chain. One can see the runaway time indicated by the vertical lines. However, the maximal discarded weight over all simulations, sites and times up to 40h/t H we measured is yet smaller than 10 −6 . Our simulations reveal that when taking a maximal Schmidt dimension of χ c = 600 instead of 300, the runaway time is shifted only about 20% to larger times and our quantities of interest are only slightly changed. Moreover, the stationary current computed with our method, which is described later, is only changed by less than 0.5%. Figure 5 shows the very good agreement between our TEBD simulations and the exact calculation using the single particle reduced density matrix (13) for spinless fermions in the tight-binding model. We use a time step δt = 0.01 and a typical simulation runs approximately 1/2 to 8 hours for a ground state calculation and about 4 hours to (seldomly) 1 week for a real time evolution. We use a multi-threaded version of TEBD with an extremely low overhead (less than 1% for 10 processors and less than 5% for up to 96 processors). The overall cost of a TEBD simulation is slightly higher than that of a DMRG calculation on a single processor. If several processors are used in parallel the overhead of DMRG calculations become rapidly prohibitive, exceeding 100 % for as few as 4 processors. 33 Thus TEBD is already twice as fast as DMRG for 4 processors and the difference is likely to increase rapidly with the number of processors. We have used 12 threads for a TEBD simulation on average.

III. CLASSICAL LC LINE
In the following, our attention will be put on a classical model: the so-called LC line, shown in Fig. 6. The LC line has been under research as an electrical circuit implementation of a Toda chain 40-42 , a nonlinear oscillator chain that, among others, describes the propagation of soliton waves. Further investigations of the LC line have so far been focused on the continuous and linear case. 43 While this setup has therefore been analyzed for the nonlinear and the linear and continuous case, our focus lies on a discrete and linear model for finite and infinite sizes.
In the following, we show that the LC line serves well as the classical representation of a one-dimensional quantum wire described by the tight-binding Hamiltonian, while our main focus lies on the behaviour of the current, especially the influence of finite-size effects and stationary values for the infinite-size limit. We present solutions for a setup in which an initial imbalance of charge carriers on the condensators leads to an oscillating, rectangular current curve as the charge carriers move back and forth in the LC line: the formerly described setup (I). Scaling the system size N to infinity leads to an enlargement of the square wave period, and for large times a stationary currentĪ is flowing. Subsequently, a method is presented to compute the stationary currentĪ for infinitely large systems using finite-size simulation results. Based on the comparison of the classical and quantum mechanical model this approach will be applied to quantum systems. Figure 6 shows the classical, electronic setup of a onedimensional wire that possesses many properties of a onedimensional quantum chain. This LC line is a combination of condensators and inductors whereas the present model is extended to enable on-site applied potentials φ i . For this setup, the following relations can be derived using elementary electrotechnical relations and Kirchhoff's laws (29) where I 0 = I N = 0 was defined. V i in equation (29) results from the potential difference between the left (φ i +U i ) and right (φ i+1 +U i+1 ) connection of the i-th inductor. C i denotes the capacity of the i-th capacitor, U i the voltage drop over the capacitor, L i the inductance of the i-th inductor and I i and J i are the currents as shown in Fig. 6. Q i is the charge of the i-th capacitor and the φ i denote externally applied potentials. From the equations above follows¨ with We consider an LC line with equal capacities and inductances and temporally constant applied potentials˙ Φ = 0 for which it follows and With Q 0,i := Q i (t = 0) and φ 0,i := φ i (t = 0) one haṡ and thus A. Stationary current In the following, we assume that initially (at t = 0) no current is flowing and that the chain is divided into two halves in which the charge carriers are uniformly distributed, i.e.
which corresponds to the previously mentioned setup (I). Figure 7 shows the current through the inductor L N/2 according to the general formula for the current at site N/2 for any even N with All values for the classical system (Q i , t, L, C, etc.) are given in S.I. units unless otherwise specified. While the rectangular shape of the dominant oscillation is understandable through the picture of charge carriers moving back and forth in the line, the rapid oscillation on top of the square wave is harder to explain and we will later take a closer look at it. The period T max cl of the square wave can be calculated by using equation (39) and (40) and taking only the smallest frequency into account (41) For large system sizes (N ≫ 1) this period is approximately given by and is thus linearly dependent on N . In the second half period of the current in Fig. 7 one can see a beat upon the rapid oscillation. This is due to the unequally distributed charges which form an oscillating pattern and so create additional frequencies in the current. This effect also occurs in quantum wires, not only after the reflection at the right border but already from the beginning due to using an initial (ground) state with non-uniform local densities.
In contrast to the oscillating current in a finite system where the charge carriers are scattered at the borders, we assume that the current in an infinitely large system is constantly flowing in one direction. To validate this assumption one has to execute the limit N → ∞ first, and then take the limit t → ∞. In the limit N → ∞ one gets from equation (39) with a = 2t/ √ LC and ϕ = 2/ √ LC. The solution is given by where J n (x) are the Bessel functions of the first kind and H n (x) denote the Struve functions. 44 The current for the infinite system (44) is shown in Fig. 8 and it shows that for smaller times the curve for finite system sizes matches the one for the infinitely large LC line. Hence, a very accurate approximation of the value of a stationary current in an infinitely large system can be extracted from the value of a corresponding finite system result. This behaviour will be later used for the analysis of quantum systems. Using the approximations (A1) for the Bessel and Struve functions we find the value of the stationary current in the infinitely large system The conductance G cl for the stationary case is given by with a local voltage drop V = (Q L − Q R ) /C in the middle of the chain.

B. Finite-time and finite-size effects
In order to describe finite-time effects in the infinite system we seek for a simpler current expression which does not only give the correct stationary value but also a good approximation for the short-time behaviour. Using the asymptotic series expansions (A2), I ∞ (t) can be transformed into the following form for t ≫ 1 of H 1 (x) in (A2) has an error of smaller order which presumably holds according to Ref. 45. As shown in Fig. 9, expression (47) is a pursuasive approximation even for smaller times and thus a good explanatory basis for short-time effects in the LC line. It states a general O(1/ √ t) decay of the amplitudes of the rapid oscillation in the current curve which is the maximal order of the error for an approximation of the stationary current (45) using a finite-time (i.e. a finitesize) result.
Instead of looking at the current through a single inductor, one can consider the current through two or more neighbouring ones. We will show that the short time behaviour of this current has some advantages compared to the simple I N/2 (t) curve. We define the quantity of interest for even N as I m,N/2 (t) = 1 2 I N/2 (t) + 1 4 I N/2+1 (t) + I N/2−1 (t) (48) which is equal to the expression 1 2 (I N/2 (t) + I N/2+1 (t)), whereas for the latter expression the subsequently explained effect also occurs for odd N for the substitution N/2 → (N −1)/2. In the limit N → ∞, t → ∞, I m,N/2 (t) converges against the same value as I N/2 (t) which implies that the stationary current is as well assessable via the current through two or more adjacent inductors. From (33) one gets with expression (36) for b k and the step distribution for the charges (38). In the limit N → ∞ it follows The integral can be solved in analogy to (43) with the substitution ξ = sin (x) and gives This result, plotted in Fig. 10 together with equation (44) for the current I ∞ (t) through a single inductor, shows almost no rapid oscillations. An approximative expression for the rapid oscillations I r (t) ≈ I dev (t) can therefore be derived using the difference between the two current expressions (51) and (44) (52) and using the approximations (A2) for the remaining Bessel functions one has Hence the rapid oscillations have a radial frequency ϕ = 2/ √ LC = 2π/T min cl and thus C. Determining the stationary current from finite system values There are various methods to calculate the stationary current from simulation results stemming from systems with finite size. The most obvious and easiest way to do so is to use the curve of a current through the middle inductor (or the mean value of currents through two or more adjacent inductors) and to compute the mean value over a certain time interval in the quasi-stationary regime. In contrast to this approach, we will present a more analytical approach using the Fourier transformation. Although one cannot prove that the following strategy can be applied to quantum mechanical systems with complicated interactions, several numerical and theoretical indications are given in sections IV and V. We start with a square wave described by its Fourier series where m determines the number of harmonics, A denotes the value of the quasi-stationary plateau and χ = 2π/T max R where T max R is the period of the square wave. This signal has certain similarities with the current curve I N/2 (t). It possesses a quasi-stationary regime which becomes infinite for χ → 0 (equals N → ∞ for the LC line) and shows a rapid and decaying oscillation on top of it. The main difference between the composition of the current curve (39) and the square wave is the presence of sine-functions that enclose the (2k − 1)-terms.
The current in (39) is dependent on η k and it is thus convenient to write I ). Choosing an η k without the enclosing sine-functions gives with ϕ ′ = ϕπ/(2N ) = π/(N √ LC). Since expression (56) now describes a square wave, one can compare it with (55) and get its quasi-stationary value which is equal to the value that was derived for the stationary current in equation (45). Remarkably the enclosing sine-functions do not change the stationary value, i.e. the expression (56) and the current in an LC line have the same stationary value for N → ∞. In case of expression (56) this value is equal to the quasi-stationary value (57) which leads to the idea to calculate the quasistationary value from a given current signal in a finite system (using the procedure described below) and use it as approximation for the stationary current value for N → ∞.
Applying the Fourier-transformatioñ term by term to expression (39) gives (40) for η k and ϕ = 2/ √ LC. If numerical data with m values are given, one can only apply a discrete Fourier transformation which preserves the area underneath a delta peak, such that for any interval ∆t containing a single delta function a · δ(x − x 0 ) the peak is of height a/∆t. Taking only the first (k = 1) and only the positive delta-peak into account and using the discrete transformatioñ and j = 0, 1, ..., m − 1, one gets for the current where Im[Ĩ N/2 (ω 1 )] is the height of the first discrete peak (k = 1). In fact, one cannot expect to know the precise form of η k for complicated quantum systems. Instead, one can take a first order approximative expression. For N ≫ 1 it holds sin π 2N ≈ π 2N , which implies that the given curve is a square wave, and thus with a relative error 44 It is important to remark that equation (62) is not only an exact result for a square wave, but also an excellent approximation for the LC line current with maximal error according to equation (63). As our simulations reveal, the charge carriers in the LC line coming from the left half first reach the right border at t ≈ T max cl /4. After this time the finite system size has a significant influence on the simulation results. Applying a discrete Fourier transformation means to assume that a periodic function is given. For that, and to minimize the upcoming error for times larger than T max cl /4 the following procedure is proposed: From a given current curve I N/2 (t) stemming from simulation results, first determine a position t = t s in the vicinity of T max cl /4but smaller than T max cl /4 -where the curve could be mirrored along a parallel to the y-axis such that the loss of continuity is minimized, which is the case for instance at inflection points. Then consider a 4t s -periodic function, constructed as follows We then apply the Fourier transformation to the constructed signal (64) and use equation (62) to calculate the corresponding stationary value of the current, see Fig. 11.

D. A final note on setups
Instead of initially dividing the system into two halves with different charges Q L and Q R , we can distribute the charge carriers homogeneously over the chain, apply a temporally constant voltage (φ 0,L − φ 0,R ) and get the same result for the current as a function of time, which can be seen in equation (36) for b k . The reason for this lies in the linear character of our model which is as a consequence independent from the chosen setup. This is the main reason why we cannot fully compare the LC line with a quantum system, since preparing the classical system in setup (II) gives the same outcomes as for setup (I), in contrast to our quantum simulation results.  Between t a and t b we observe an approximately rectangular wave with a period T max which diverges with increasing system length (the corresponding leading frequency ω max = 2π/T max converges towards 0 in the thermodynamic limit), see Fig. 12. The period of the rectangular oscillation is given by (Appendix B) For spinless fermions, using a semi-classical picture, particles first flow from one side of the system to the opposite side because of the inhomogeneous density (first setup) or the potential difference (second setup). Then they are reflected by the hard wall represented by the chain edge. As there is no dissipation in our model, all reflected particles flow back with the same velocity in the opposite direction until they reach the other chain edge and are again reflected and so on.
The progressive degradation of the rectangular signal can be understood using the same picture. First, all particles flow in the same direction but, as they have different velocities, they progressively come out of phase. For long times t ≫ t b , which can be checked up to the numerical double limit t ≈ 10 300 for non-interacting fermions (V H = 0), our simulations show that the current does not go to zero but continues to oscillate with a period T max . This can be understood in the picture of the classical LC line as well as for spinless fermions through the dominant amplitude of the ω max oscillation. It is basically the same effect as if all contributions to a series representation of a square wave (55) would 'randomly' come out of phase. In that case, the contribution with the lowest frequency ω max , which is the frequency of the rectangular oscillation, determines the frequency of the remaining irregularly shaped curve.
B. Rapid oscillation on top of the square wave Setting the two expressions (41) and (65) for the periods of the classical and quantum mechanical system equal one gets Further comparison with equation (54) gives for the period of the rapid oscillation Instead of preparing the system in a ground state of a Hamiltonian with applied potentials, one can initially decouple all sites and fill the left half of the system with (uncorrelated) particles, leaving the right half empty. This almost corresponds to choosing ∆ǫ = 4t H , insofar as a ground state with this applied voltage and coupled reservoirs still has some correlations and unequally distributed on-site particle densities in it. The resulting current curve in Fig. 12 in the first half period closely resembles the curve shown in Fig. 7 for the LC line.
The magnification in Fig. 12 shows a rapid oscillation with a period according to equation (67), in contrast to the oscillation shown in Fig. 13 for setup (I) which is 'disturbed' by oscillations stemming from correlations and an unequal distribution of particles in the ground state due to a lower applied voltage ∆ǫ < 4t H . Although the oscillation in Fig. 13 for setup (II) is very regular, it does not stem from the same 'classical' origin, which is confirmed by the fact that on the one hand two adjacent currents do not cancel out and on the other hand the oscillation does not fulfill equation (67).  Since in the following we are more interested in the stationary current flowing from an infinitely large source (the left half of the chain) to a an infinitely large drain (the other half), we reference to Ref. 46 for a discussion of time-dependent currents and to Ref. 35 for a detailed analysis of short-time and finite-size effects (in the IRLM).

C. Applicability of the extrapolation approach to quantum systems
Although we have seen the close connection between the tight-binding system for initially uncorrelated particles and the classical LC line, there exists a major difference in the behaviour of the respective currents: the damping of the rapid oscillations. While for the classical LC line a O(t − 1 2 ) dependence of the amplitudes according to equation (47) is predicted, this decrement is different for a quantum system.
The expectation value of the current for setup (I) and for initial conditions like in Fig. 12 (one completely filled and one completely empty half) is for N → ∞ given by 47,48 with ω = 2t H /h. Using the asymptotic series expansions (A2) yields for the current which is plotted together with expression (68) in Fig. 14.
The curve highly coincides with the one shown in Fig. 12 for t < T max /4 which confirms the assumption that the borders in the finite quantum system for setup (1) only significantly change the current after a time T max /4, determined by the velocity of the particle density front wave 48 which moves from the middle of the chain to the borders. The quasi-stationary valuē J = 4et H /h from (69) suits equation (10) for ∆ǫ = 4t H and the amplitudes of the rapid oscillation are damped with O(1/t). Moreover, expression (67) for the period of the rapid oscillation -which was derived only by comparison with the classical model -is confirmed by equation (69).
Since the latter analysis shows a O(1/t) damping of the rapid oscillation, the open question at hand is, if it is reasonable to use our extrapolation method for a quantum system. In fact, a pure square wave described by its Fourier series (55) has also a damping of the rapid oscillation of O(1/t) and highly resembles the current described by (68) or depicted in Fig. 12 when choosing an appropriate amplitude and number of harmonics. Since the mathematical error (63) of our method stems from the difference of the analyzed curve to a square wave, the extrapolation approach provides an even smaller error than (63) for a quantum system.

V. SPINLESS FERMION MODEL
In this section we discuss our results for the stationary current and the finite-system period in the spinless fermion model. The non-linear dynamics of this model and the closely related spin 1 2 XXZ chain have been studied in several works using td-DMRG methods 23,39,49 . Nevertheless, the steady-state transport properties of the spinless fermion model have not been investigated yet.

A. Current-voltage characteristics
Firstly, we remark that the value of the potential bias is given by the parameter ∆ǫ of the Hamiltonian from equation (4), and it is not the result of a quantum measurement of an observable. Thus it corresponds to an external field applied to the system, and in an interacting system it is not always equal to the actual potential difference which is measured in experiments. 6,50 Moreover, in experiments the sample (quantum wire) must be connected by leads to the voltage source and measurement apparatus. This modifies the low-energy, low-temperature behaviour of the sample below some crossover energy scale. 7,51 Concerning this matter, the linear conductance of a TLL attached to one-dimensional leads has been calculated in several works. [51][52][53][54] It has been shown that the conductance of this setup is given by e 2 /h (which corresponds to equation (10) with K = 1), independently of the Coulomb interaction V H within the wire, and that this result coincides with experimental outcomes.
Our analysis of on-site particle densities shows that the densities for the first and last site of the chain only change after T max /4. This is a major reason why we only take values of the current for times smaller than T max /4 into account, as described in our extrapolation approach.
To test the overall procedure for computing the stationary current we have first calculated the current-voltage characteristic of a non-interacting chain (V H = 0). This characteristic is known exactly for the first setup 47,48 : For the second setup, exact results were calculated numerically using the one-particle equation of motion with a system size of N = 1000. We have found that our procedure yields stationary current values which agree with an overall error of less than 5% with the exact results. As a second test we have calculated the current-voltage characteristic of the IRLM. The generic shape of this curve is known from a field-theoretical analysis, and highly accurate td-DMRG results are available for a quantitative comparison. 12 In that regard, we have found that the results from our method agree very well with the numerical outcomes from Ref. 12. We nevertheless have to remark that we found the IRLM far easier to simulate than the spinless fermion model, since for higher |V H | the correlations (and thus the needed Schmidt dimension) within the spinless fermion model grow more quickly than in the IRLM. As already mentioned in section II D the parameters used for our simulations can result in larger errors for high voltages and Coulomb interactions in the TEBD data for the current as a function of time. Additionally, finitesize effects also grow with |V H | and ∆ǫ. This can lead to less accurate extrapolations for the steady-state current in that regime than for non-interacting fermions and the IRLM.
We now discuss the current-voltage characteristics which have been obtained with the methods described above. In setup (I) the current of the non-interacting system (V H = 0) is strictly proportional to ∆ǫ = |eV | up to the band width 4t H and then remains constant. 47 For V H > 0 we see in Fig. 15 that the current increases sub linearly with ∆ǫ for a fixed interaction strength V H and that it decreases monotonically with increasing V H for a fixed potential bias ∆ǫ. The magnification in Fig. 16 shows a comparison with the linear response in the TLL theory (10) for small ∆ǫ, using the Bethe Ansatz solution parameters (3). The good agreement confirms the validity of our approach. Obviously, an increasing V H does not only reduce the current but also the range of the potential bias ∆ǫ for which the linear response approximation is accurate. In Fig. 17 we observe a behaviour for attractive interactions V H < 0 which is similar to the non-interacting case. The current increases almost linearly with ∆ǫ up to some V H -dependent threshold potential where it saturates at its maximal value. Increasing |V H | scales the maximal current down but causes a higher linear conductance for small ∆ǫ. Figure 18 shows that for negative interactions there is also an excellent agreement between our results and the linear response in the TLL theory (10) for small ∆ǫ. Again increasing |V H | seems to reduce the range of the potential bias ∆ǫ for which the linear response approximation is accurate. The occurrence of current plateaus for large potential differences can be easily understood. In setup (I) one half of the chain initially contains more particles than the other one due to the applied potential bias ∆ǫ. When this bias is large enough the initial state consists of one completely filled and one completely empty half and it does not change for higher ∆ǫ. In the non-interacting case this saturation occurs exactly at ∆ǫ = 4t H . 47,48 When we add an attractive interaction V H < 0 the particles are more likely to stick together and therefore even more particles gather in one half for the same applied potential. Thus, saturation (i.e., an initial state with one completely filled and one completely empty half) is reached for smaller ∆ǫ as |V H | increases as seen in Fig. 17. We do not understand yet why the plateau heights (i.e. the maximal current) are lowered by increasing |V H |. For a repulsive interaction V H the effect is reversed and less particles gather in one half of the chain for a given potential difference when V H increases. Thus we expect that saturation occurs at higher values ∆ǫ > 4t H beyond the potential range shown in Fig. 15. We can easily determine for which potential difference ∆ǫ saturation occurs. Removing a single particle from the completely filled reservoir or adding one particle to the empty half costs an energy ∆ε/2 − V H − 2t H . Thus saturation occurs if (70) Figure 17 shows that this approximation fits well to the numercial data for V H ≤ 0 whereas the saturated regime according to (70) for V H > 0 lies outside the potential range shown in Fig. 15, as mentioned above. While for setup (I) the differential conductance is always positive (or zero after saturation), in setup (II) we observe a negative differential conductance. Figure 19 shows our results for setup (II) with V H ≥ 0. The current seems to vanish for very large potential biases as predicted by strong-coupling perturbation theory for this setup. We note that the current becomes systematically weaker with increasing V H for a fixed small potential bias. For larger ∆ǫ the behaviour of I as a function of V H is no longer monotonic. In addition, we see that the onset of the negative differential conductance (i.e., the position of the maximum of I as a function of ∆ǫ) shifts to higher values with increasing V H . The magnification for small ∆ǫ, Fig. 20, confirms again that our results agree with the TLL theory in the linear regime and that the range of the linear response regime shrinks with increasing interaction strength.
Results for attractive interactions V H < 0 are shown in Fig. 21. There we clearly observe that the current vanishes for large potential differences as predicted by strong-coupling perturbation theory for setup (II). We also see that the current becomes systematically larger with increasing interaction strength |V H | for a fixed small potential bias. For intermediate ∆ǫ the behaviour of I as a function of V H is not monotonic and for large enough ∆ǫ > ∼ 2t H the current decreases systematically with in-  creasing interaction strength. Moreover, we see that the onset of the negative differential conductance occurs at lower values of ∆ǫ for increasing |V H |. Finally, Fig. 22 shows that in this case our results also agree with the TLL theory in the linear regime, but it seems that the interaction reduces rapidly the range of potential biases for which the linear response is accurate. The negative differential conductance which is observed in setup (II) for all V H can be understood in terms of the overlap of the energy bands for elementary excitations. Our simulations show that there is no significant difference in the current curve whether one initially does or does not couple source and drain (i.e., the two halves of the system). For the following considerations we can therefore regard source and drain as initially decoupled, having two separate but equal energy bands. Then we apply a step-like voltage in setup (II) and thus shift the two bands against each other. For V H = 0 only overlapping occupied one-particle states in the source and unoccupied ones in the drain contribute to the current. While this overlap rises with increasing |∆ǫ| from 0 to 2t H , it diminishes from |∆ǫ| = 2t H and reaches zero at 4t H . As a result, the current is maximal for ∆ǫ ≈ 2t H and vanishes for large ∆ǫ. Similarly, we can understand the non-monotonic behaviour in interacting cases (V H = 0) in terms of the overlap of the elementary excitation bands in the spinless Luttinger liquids in the two halves of the system. However, the effective bandwidth is renormalized exactly as the Fermi velocity in equation (2). Therefore, the maximum of the current is reached for a smaller potential bias ∆ǫ when V H becomes negative as shown in Fig. 21, and shifted to a higher potential bias when V H increases as seen in Fig. 19. Our conclusion agrees with the findings in Ref. 36 where it is shown that within a similar one-dimensional model a negative differential conductance is mainly caused by finite electrode bandwidths.
B. Influence of VH on the period T max in the finite quantum system Our results show a further effect of V H on the finitesystem current. While for V H > 0 the period of the rectangular oscillation becomes smaller, it grows for negative interactions. Figure 23 shows the period T max of the rectangular current oscillation in the finite system as a function of the inverse renormalized Fermi velocity v −1 from expression (2). The T max values are mean values from various ∆ǫ for each setup. These values agree perfectly with the dashed lines in Fig. 23 which are given by the equation Thus, the period is fully determined by the time-scale of the non-interacting system T max (V H = 0) from (B4) and the renormalized Fermi velocity v from expression (2).

VI. SUMMARY AND CONCLUSION
We have presented a method for investigating the non-linear steady-state transport properties of onedimensional correlated quantum systems. First, we have analyzed finite-size effects and transient currents in a related classical model, the so-called LC line. This exact analysis shows how steady-state currents in the thermodynamic limit can be calculated from transient currents in finite systems. We have found that for currents in quantum systems, in particular the spinless fermion model, finite-size effects can be understood from the behaviour of the LC line.
On the basis of the strong connection between the finite-size behaviour in classical and quantum models we have developed a method to extract a stationary current from simulations of finite-size quantum systems. We only need numerical data I(t) up to a time of the order of T max /4 where T max ∝ N is the period of the approximately rectangular signal I(t) in a quantum system of size N . The non-equilibrium dynamics of correlated quantum systems has been calculated using the timeevolving block decimation method (TEBD). We have implemented a multi-threaded version of TEBD with an extremely low overhead (less than 1% for 10 threads) which allows us to simulate the non-equilibrium dynamics up to the time-scale 2T max .
We have determined the full I-V characteristic of the spinless fermion model with nearest-neighbour hopping t H and interaction V H using two different setups to generate currents. In setup (I) the initial state has different particle numbers in its two halves due to an applied potential difference while the system evolves in time with an overall equal on-site potential. In setup (II) we calculate the initial state without a potential difference but turn it on for the real time evolution. For non-interacting fermions (V H = 0) our outcomes agree with the analytical solutions and with the results of the one-particle equation of motion method. Additionally, we have checked that our results coincide with the field-theoretical analysis and td-DMRG simulations 12 for the IRLM.
For interacting fermions we have found that the steadystate current is independent from the setup in the linear regime |eV | ≪ 4t H . Moreover, our numerical data agree with the predictions of the Luttinger liquid theory combined with the Bethe Ansatz solution. For larger potentials V the steady-state current depends on the currentgenerating setup. This difference is well understood for non-interacting systems 35,36 but can only be explained qualitatively for our interacting system. With setup (I) we have found that the current increases with V and saturates at a finite value when the potential difference is high enough to separate the initial state in one filled and one empty half. With setup (II) we observe a negative differential conductance that can be understood in terms of the overlap of the elementary excitation bands in the spinless Luttinger liquids in the two halves of the system. Both effects -the plateaus and the negative differential conductance at large potential bias -are due to the finite bandwidth of the system which agrees with the findings in Ref. 36. However, in our case the interaction renormalizes the effective bandwidth exactly as the Fermi velocity in equation (2) so that the current maxima and the onset of the plateaus depend on the strength of the interaction V H . The change in the typical time-scale T max also solely depends on the renormalized Fermi velocity.
The methods presented in this work can be applied to other systems, such as models of quantum dots or wires coupled to leads including electronic and also bosonic degrees of freedom. We have already checked the applicability of these methods to the IRLM, the one-dimensional Hubbard model away from half-filling and also two-leg ladder systems. Hence, we believe that they will be very useful to study non-linear transport properties of correlated low-dimensional conductors. Although we have mostly tested our appraoch in the limit of transparent coupling between source and drain, we think that it can also be applied when the hybridization between leads is very small. This issue will be examined in future works.