Landauer transport as a quasisteady state on finite chains under unitary quantum dynamics
J. P. Santos Pires B. Amorim J. M. Viana Parente Lopes

#1#1

#1|#1|

#1⟨#1⟩

|#⟩1|#1⟩

⟨#|1⟨#1|

⊗

⟨#|1⟩#2⟨#1|#2⟩

Ω

inf∞



up201201453@fc.up.pt

Centro de Fsica das Universidades do Minho e Porto 

Departamento de Fsica e Astronomia, Faculdade de Cincias, Universidade
do Porto, 4169-007 Porto, Portugal

amorim.bac@gmail.com

Centro de Fsica das Universidades do Minho e Porto 

University of Minho, Campus of Gualtar, 4710-057, Braga, Portugal

jlopes@fc.up.pt

Centro de Fsica das Universidades do Minho e Porto 

Departamento de Fsica e Astronomia, Faculdade de Cincias, Universidade
do Porto, 4169-007 Porto, Portugal



(Accepted
in the Physical Review B - March 2020)



In this paper, we study the emergence of a Landauer transport regime
from the quantum-mechanical dynamics of free electrons in a disordered
tight-binding chain, which is coupled to finite leads with open boundaries.
Both partitioned and partition-free
initial conditions are analyzed and seen to give rise, for large enough
leads, to the same spatially uniform quasi-steady-state current, which
agrees with the Landauer value. The quasi-steady-state regime is preceded
by a transient regime, which lasts for a time proportional to the
length of the disordered sample, and followed by recursions, after
a time that is proportional to the lead size. These theoretical predictions
may be of interest to future experiments on transport of fermionic
ultra-cold atoms across optical lattices. We also observe finite-size
current oscillations, superimposed on the quasi-steady state, whose
behavior depends crucially on the conditions initially imposed on
the system. Finally, we show how a time-resolved Kubo formula is able
to reproduce this Landauer transport regime, as the leads grow bigger.



Landauer transport as a quasisteady state on finite chains under unitary
quantum dynamics
    J. M. Viana Parente Lopes
    December 30, 2023
=========================================================================================





§ INTRODUCTION


The study of electronic transport is amongst the main goals of condensed
matter physics. In the regime of small length scales and low temperatures,
the mesoscopic transport regime, quantum coherence effects play a
dominant role in the propagation of electron states. In such a case,
transport can no longer be seen as a bulk phenomenon, but instead
depends on device-specific details such as the geometry of and nature
of the electrodes, as well as the specific distribution of disorder
in the sample.

A theoretical description of mesoscopic transport was first developed
by Landauer <cit.> and later generalized
by Bttiker <cit.>. In the now called
Landauer-Bttiker formalism, the problem of stationary mesoscopic
transport is recast as a scattering problem, where single-electron
states incoming from the leads are transmitted across a central device.
The current may then be expressed as a sum over the transmission probabilities
of the occupied incoming lead states. In parallel to this work, Caroli
<cit.> applied the non-equilibrium Green's function
formalism of Kadanoff-Baym <cit.> and Keldysh
<cit.> to the calculation of mesoscopic transport.
The obtained expression has a structure similar to the Landauer-Bttiker
(LB) one, but with the transmission coefficient now expressed in terms
of Green's functions of the central device and spectral functions
of the leads. While apparently distinct, the two approaches lead to
the same result, as implied by the Fisher-Lee relation <cit.>
between transmission coefficients and Green's functions (for a detailed
proof, see Wimmer <cit.>). Central to both approaches
are the assumptions that the leads attached to the central device
are semi-infinite and the occupation of the incoming single-electron
states is determined by independent Fermi energies on each lead. Moreover,
both methods are only able to describe steady-state transport. We
also point out that the Landauer formula has also been derived within
the theory of non-equilibrium steady states <cit.>.

If one is interested in the transient dynamics and how this steady
state is reached, the matter of what the initial condition of the
system was, becomes relevant. At the theoretical level, two initial
conditions have been historically considered: (i) In the
partitioned approach <cit.>,
the leads and the central device are assumed to be initially disconnected,
each being in equilibrium with independent Fermi levels. This Fermi-level
imbalance takes into account the bias applied to the mesoscopic device.
Then, the leads and the device are suddenly brought into contact allowing
a charge current to flow. (ii) In the partition-free
approach <cit.>,
the leads and the device are assumed to be connected from the beginning
and in global equilibrium with a common Fermi energy. Then, a potential
bias between the leads is suddenly applied to the connected system.
It has been shown, for the case of a single-level central device,
that the same steady-state current is reached in both approaches,
for the same time-dependent perturbation <cit.>.
Furthermore, the value of the steady-state current does not depend
on the history of the time-dependent perturbation, provided it reaches
the same constant value in the future. Crucial to this result is the
fact that the leads have a continuum spectrum (as it occurs for semi-infinite
leads), which allows a loss of memory about the initial conditions,
provided there are no bound states in the central device. The existence
of bound states inside the device is known to cause persistent current
oscillations <cit.>.
We also point out that, under certain circumstances, interactions
might prevent the formation of a steady state <cit.>.
Finally, the equivalence of the steady-state current reached in the
partitioned and partition-free cases was shown,
with rather broad assumptions and including the presence of interactions,
in a mathematically rigorous way in Refs.<cit.>.

Very recently, Purkayastha et
al <cit.> also
discussed the problem of relating the linear response theory in a
general interacting system coupled infinite leads to the open system's
Kubo formula. The authors looked at the asymptotic behavior of the
integrated current-current correlators at different times, in the
limit of very large systems and times. This way, they found a surprising
disagreement between the transport classification given by both approaches
<cit.> for the critical one-dimensional Aubry-Andr-Harper
model. This was attributed to the non-commutativity of the two referred
limits. Nevertheless, their analysis was always done in a partitioned
setup and in linear response theory.

In more recent years, a significant effort was devoted to the study
of time-dependent transport and transient dynamics in mesoscopic systems
attached to infinite leads <cit.>.
However, the investigation of time-dependent transport in systems
where the leads are finite (but possibly very large) has received
much less attention. Initial work on this problem was made in Refs.
<cit.>,
in which a micro-canonical method was developed to deal with quasi-steady-state
transport in finite systems, where the leads are initially connected
to the device, but one of them is partially or fully depleted of particles.
More recently, the same problem was considered in Pal et al
<cit.>, which studied time-dependent transport
through a quantum dot connected to two systems with a quasi-continuum
spectrum (discrete, but dense), which take the role of finite leads,
in the partitioned approach. In all of these works, it was
observed that after the transients died out, a quasisteady state transport
regime emerges. However, the dependence of the quasisteady state on
the system initial condition has never been explored so far. Is the
quasi-steady-state current independent of the initial preparation
of the system? Does the manner in which the quasisteady state is reached
depend on these initial conditions? These are the questions we try
to answer in this work. Notice that these issues are not of mere theoretical
interest. It is true that for realistic electronic devices, consideration
of finite leads is somewhat artificial.[Unless one is interested in times comparable to the discharge time
of a battery, the electrodes and voltage source that are connected
to a mesoscopic device can be seen as infinite.] However, cold atoms trapped in optical lattices have emerged as a
platform for the experimental study of transport properties <cit.>
[also studied theoretically in Chien et al. <cit.>].
In these systems, due to the limited sizes of optical lattices, the
“leads” are necessarily finite. Furthermore, while in realistic
low-dimensional electronic devices, one expects electron-electron
interactions play a significant role <cit.>,
in optical lattices the inter-particle interaction can be tuned down
to zero, thus enabling a proper study of transport for non-interacting
fermions <cit.>. As such, we
expect that our theoretical results will be of experimental relevance
for transport in ultra-cold-atoms setups.

The purpose of this work is to further explore how a steady-state
transport regime emerges from quantum time-evolution in non-interacting,
fermionic systems with finite, but large leads, and how the initial
state of the system affects this process. By combining numerical and
analytical work, we study the time-dependent current dynamics in a
prototypical one-dimensional non-interacting tight-binding model with
disorder, analyzing in detail how the current dynamics depends on
the initial conditions (partitioned vs partition-free)
and on the size of the finite leads. We employ a full quantum time
evolution, starting from both initial conditions, to study the time-dependent
current upon the sudden connection of the appropriate perturbation.
For the partition-free case, we also derive a time-dependent
Kubo formula, which allows us to see rigorously how a linear Landauer-Bttiker
formula, involving only quantum transmittances, emerges from an unitary
time evolution in the limit of very large leads.

The text is organized as follows. In Sec. <ref>,
we introduce the one-dimensional tight-binding model Hamiltonian that
will be used throughout the rest of the paper and detail both the
partitioned and partition-free approaches. In Sec. <ref>,
we describe the numerical methods used for calculating the time-dependent
local current from the unitary dynamics of the finite system and also
the steady-state Landauer current for infinite leads. The main numerical
results are then presented in Sec. <ref>, where
the time-evolution of the non-equilibrium current is systematically
analyzed as a function of the bias, the size of the finite leads,
and the central sample's disorder and size. Finally, in Sec. <ref>,
we provide analytical insight into the numerical results of Sec. <ref>,
by developing a time-dependent Kubo formula for the partition-free
approach and expressing it in terms of complex reflection and transmission
coefficients of the central sample. In Sec. <ref>,
we discuss the obtained results and conclude the paper.






§ MODEL HAMILTONIAN AND INITIAL CONDITIONS


We will study the current dynamics of non-interacting electrons in
a finite one-dimensional tight-binding model, with nearest-neighbor
hoppings. The tight-binding chain is composed by a total of L sites,
with the central L_s sites, the sample, having an on-site Anderson
disorder and being subject to a constant electric field. The sites
outside the sample region form the left and right leads [each with
L_l=(L-L_s)/2 sites], are not disordered and hold
a constant electrostatic potential. They will refer to the different
regions in the chain as left lead (LL), sample (S), and right lead
(RL). An illustrative scheme of this setup is shown in Fig. <ref>.
For times t>0, the dynamics of the system is governed by the time-independent
Hamiltonian






    ℋ(t>0)=∑_n=0^L-1(ϵ_n^d-ev_n^e)c_n^†c_n-w∑_n=0^L-2(c_n+1^†c_n+c_n^†c_n+1),

where c_n^†(c_n) are creation (annihilation)
operators for an electron at the chain site n, w is the nearest-neighbor
hopping amplitude, e>0 is the fundamental charge and v_n^e
is the electrostatic potential. According to the previous dis-









-cussion v_n^e has the form

    v_n^e=Δ V/2    , n∈0,...,L_l-1
    (1/2-n-L_l+1/L_s+1)Δ V     , L_l≤ n<L_l+L_s
    
    -Δ V/2    , n∈ L_l+L_s,...,L
    ,

where Δ V is the applied potential bias, and ϵ_n^d
is the Anderson on-site potential disorder, which is only present
in the sample sites, and are taken as random numbers uniformly distributed
inside [-W/2,W/2].

We will study the current dynamics in this system both in the partitioned
and partition-free approaches. In both cases, the dynamics
for t>0 are governed by the Hamiltonian of Eq. (<ref>),
with only the initial state being different.

In the partitioned approach, the initial state is formed
by occupied states for the partitioned system, with the bias already
applied. The partitioned system is described by the Hamiltonian,





    ℋ^P(t=0)=ℋ_LL^P+ℋ_S^P+ℋ_RL^P,

with ℋ_LL, ℋ_S, and ℋ_RL
the Hamiltonians, respectively, for the decoupled left lead, sample,
and right lead, to wit

    ℋ_LL^P    =∑_n=0^L_l-1(-ev_n^e)c_n^†c_n-w∑_n=0^L_l-2(c_n+1^†c_n+H.c.),
    ℋ_S^P    =     ∑_n=L_l^L_l+L_s-1        (ϵ_n^d-ev_n^e)c_n^†c_n-w      ∑_n=L_l^L_l+L_s-2             (c_n+1^†c_n+H.c.),
    ℋ_RL^P    =       ∑_n=L_l+L_s^L-1    (-ev_n^e)c_n^†c_n-w      ∑_n=L_l+L_s^L-2             (c_n+1^†c_n+H.c.).

The initial occupation of the single-electron states is determined
by the independent Fermi levels for each region. Hence, we write ε_F,LL=ε_F+Δ V/2,
ε_F,S=ε_F and ε_F,RL=ε_F-Δ V/2,
as the chemical potential for the left lead, central sample, and right
lead, respectively. ε_F is a reference chemical
potential. The initial state is thus described by the reduced density
matrix




    ρ^P(t=0)=∑_r=LL,S,RL∑_α_rf_r,α_r^P|Ψ_r,α_r^P⟩⟨Ψ_r,α_r^P|,

where |Ψ_r,α_r^P⟩ are the
independent single-electron eigenstates of the initial partitioned
Hamiltonian belonging to region r, ℋ_r^P,
with an energy ε_r,α_r^P. At any temperature,
the initial occupation of the states is given by the factor f_r,α_r^P=f(ε_r,α_r^P-ε_F,r),
with r=LL,S,RL and f(ε)=(e^βε+1)^-1
being the Fermi-Dirac distribution function. Throughout this work,
we will restrict ourselves to the T=0 case, where f(ε)=Θ(-ε)
and Θ(x) being the usual Heaviside step function.
The hoppings between the leads and the sample are then suddenly switched
on and the time evolution of these states is generated by the Hamiltonian
in Eq. (<ref>).

In the partition-free approach, the contact between the sample
and the leads is already established in the initial state, but the
bias is not yet applied. Therefore, the initial condition is determined
by populating the eigenstates of the partition-free Hamiltonian

    ℋ^PF(t=0)=∑_n=0^L-1ϵ_n^dc_n^†c_n-w∑_n=0^L-2(c_n+1^†c_n+H.c.),

up to a commonly defined Fermi energy, ε_F.
This initial state is thus described by the reduced density matrix

    ρ^PF(t=0)=∑_αf_α^PF|Ψ_α^PF⟩⟨Ψ_α^PF|,

with |Ψ_α^PF⟩ being eigenstates
of Eq. (<ref>) with an eigenenergy
ε_α^PF. The initial occupation factor
of the states is similarly given by f_α^PF=f(ε_α^PF-ε_F).
In this case, the sudden perturbation driving the current is the connection
of the bias potential, v_n^e, at t=0, after which
the time evolution is governed by the Hamiltonian of Eq. (<ref>).

We end this section, by noting that the charge current flowing from
site n to site n+1, for the Hamiltonian of Eq. (<ref>),
is given by




    ℐ^n=ew/iħ(c_n+1^†c_n-c_n^†c_n+1).








§ NUMERICAL METHODS FOR CURRENT EVALUATION




 §.§ Time-resolved current from quantum
evolution of eigenstates


The dynamics of the system, in the partitioned approach after suddenly
switching on the lead-sample hoppings, or in the partition-free approach
after suddenly switching on the external bias, is governed by the
Hamiltonian (<ref>). Therefore, in
both approaches and for t>0, the reduced density matrix of the
system evolves according to





    iħdρ(t)/dt=[ℋ(t>0),ρ(t)].

The solution for this equation, with initial condition given by either
Eq. (<ref>) or (<ref>),
is given by





    ρ(t)=∑_αf_α|Ψ_α(t)⟩⟨Ψ_α(t)|,

with the single-electron states evolving according to Eq. (<ref>):
|Ψ_α(t)⟩ =e^-i/ħℋ(t>0)t|Ψ_α⟩,
with |Ψ_α⟩ being the single-electron
eigenstates of either ℋ^P(t=0) or ℋ^PF(t=0)
(with occupation f_α), for the partitioned and
partition-free approaches, respectively.

The expected value of the current, as a function of time, is given
by




    I^n(t)=ew/iħ∑_α∈occupied(⟨.Ψ_α(t)|n+1⟩⟨ n|Ψ_α(t).⟩.
    .-⟨.Ψ_α(t)|n⟩⟨ n+1|Ψ_α(t).⟩),

where |n⟩ represents the state localized at site
n. We also used the fact that at t=0, f_α=1 for initially
occupied states and f_α=0 for empty states. The above expression,
allows us to evaluate the current flowing from site n to site n+1
provided we know the time evolution of the initial single-electron
states. Although correct, Eq. (<ref>) is
not very convenient from a numerical point of view, since for each
initially occupied state, we would have to perform one time-evolution.
A more convenient expression is obtained by writing ⟨.Ψ_α(t)|n⟩ =⟨Ψ_α|e^i/ħℋ(t>0)t|n⟩ =⟨.Ψ_α|n(-t)⟩,
such that instead of evolving the initial eigenstate forwards in time,
we evolve the localize states backwards in time.[Note that, despite being an evolution for negative times, the time-evolution
operator used to do it is the one which includes the external perturbation:
lead-sample hopping in the partitioned approach, and bias in the partition-free
approach.] 

The current can therefore be expressed as


    I^n(t)=2ew/ħIm∑_α⟨ n(-t)|Ψ_α.⟩⟨.Ψ_α|n+1(-t)⟩ .


Despite being equivalent to Eq. (<ref>),
this last expression allows for a great gain in computational efficiency,
as the the number of required time evolutions is reduced from 𝒪(L)
to only two, for each single-time calculation of the current between
sites n and n+1.

Numerically, the time-evolution of the localized states |n⟩,
for very large chains, is computed efficiently using a polynomial
Chebyshev expansion <cit.>
of the time evolution operator 𝒰(t)=e^-i/ħℋt
(for details, see Appendix <ref>). Finally,
the single-electron eigenstates of ℋ(t=0)
were calculated “on the fly” using a memory-efficient algorithm
developed by Fernando <cit.>.



 §.§ Landauer-Bttiker formula for
the steady-state current


In Sec. <ref>, the time-resolved current across
the sample with finite leads, calculated using the quantum evolution
of the occupied states, will be compared with the steady-state value
for the current as given by the Landauer-Bttiker formula for the
same sample attached to infinite leads. We evaluate this current using
the Caroli-Meir-Wingreen form of the Landauer-Bttiker formula <cit.>,
which for a two-terminal device, at zero temperature, reads as 





    I_LB=e/2πħ∫_ε_F-eΔ V/2^ε_F+eΔ V/2dε T(ε),


where the energy dependent transmission function is expressed
in terms of Green's functions as





    T(ε)=Tr[G^A(ε)·Γ_RL(ε)·G^R(ε)·Γ_LL(ε)],

where Γ_LL/RL(ε) are real-space
spectral functions of the unattached leads and G^R/A(ε)
are the real-space retarded/advanced Green's function of the central
sample, in the presence of the leads. For our particular one-dimensional
model, the leads' spectral functions are matrices with the only non-zero
elements between boundary sites, i.e., Γ_LL^L_l-1,L_l-1(ε)=Γ_LL(ε)
and Γ_RL^L_l+L_s,L_l+L_s(ε)=Γ_LL(ε)
, where Γ_LL/RL(ε)=w^2ρ_LL/RL(ε).
The functions ρ_LL/RL(ε) are surface
density of states of the leads which may be computed analytically
yielding:





    ρ_LL/RL(ε)=Θ(4w^2-(ε∓eΔ V/2)^2)
    ×1/w^2√(4w^2-(ε∓eΔ V/2)^2).

Therefore, the final form for the transmission function reads


    T(ε)=w^4ρ_LL(ε)ρ_RL(ε)|G_L_l-1,L_l+L_s^R(ε)|^2.

For each central sample, the retarded Green's function in Eq. (<ref>)
was calculated by using the well-known recursive Green's function
method <cit.>,
using the surface Green function of the semi-infinite one-dimensional
leads as boundary conditions, as detailed in Appendix <ref>.



§ NUMERICAL RESULTS AND COMPARISON WITH
THE LANDAUER FORMULA


We evaluated the time-dependent current in finite open chains, using
the method described in Sec. <ref>,
for both clean and disordered samples and considering both the partitioned
and partition-free initial conditions. This current was then
compared with the Landauer expression for the steady-state current
flowing through the same sample attached to infinite leads, as described
in Sec. <ref>.

Our results are summarized in Fig. (<ref>).
For both initial conditions, three transport regimes are clearly distinguished
for large enough leads: (i) initially, we have a transient
regime up to a time t_stab, after which (ii)
the current tends to an approximately constant quasi-steady-state
value, which last up to (iii) a recurrence time, t_r,
after which an inversion of the current occurs. These three transport
regimes have been previously reported for non-interacting fermions
in systems with finite leads, in Refs. <cit.>,
with an initial state where the leads are connected to the device,
but one of the leads is depleted. These transport regimes have also
been observed for interacting fermions (at a TDDFT-ALDA level) with
finite leads by Bushong et al <cit.>.
These three regimes have also been discussed in Di Ventra et
al <cit.> and Pal et al <cit.>.
However, to the best of our knowledge the effect of different initial
conditions, in each of the transport regimes has never been explored.

In the following, we will analyze in detail each of these regimes,
analyzing how the different time-scales depend on the variables of
the problem.










 §.§ Transient Behavior and Stabilization Times


As one could expect, the transient behavior depends on the initial
preparation of the system, being different for the partitioned
and partition-free approaches, as evident in Fig. <ref>.

In the partitioned case, as can be seen in Figs. <ref>
(and also in Figs.  <ref> and <ref>),
the current is initially close to zero up to a build-up time, t_b,
after which the current dramatically increases, overshooting the Landauer
current value. We interpret t_b as the time it takes for
a fermion close to the initial reference Fermi level, ε_F,
to travel from the lead-sample boundaries to the hopping where the
current is being probed. We will define L_meas, L and
L_meas, R as the distance from the point where the current
is being measured in the sample to the left and right lead boundaries,
respectively. We will refer to the smallest of these distances as
L_meas=min(L_meas, L,L_meas, R).
According to this interpretation, the build-up time is given by t_b=L_meas/v_F
where v_F is the velocity of a state initially at the
reference Fermi level (the initial Fermi level defined for the sample,
for the partitioned case, and the initial global Fermi level in the
partition-free case). In Fig. <ref>, we show the
current for different samples, measured at L_meas, L=64.
The dashed-dotted vertical lines in Figs. <ref>
and <ref> indicate the time t_b=L_meas/v_F.
The coincidence of these lines with the 



sharp rise of the current confirms our interpretation. This
is further confirmed in the top panel of Fig. <ref>,
where the current measured at different sites is shown as a function
of time, with the dash-dotted vertical lines indicating the times
t_b^L/R=L_meas, L/R/v_F. Once
again, we can see that the current in the partitioned case is nearly
zero up to the t_b=min(t_b^L,t_b^R).

In the partition-free setup one observes a gradual increase
of the current from the beginning. As can be seen in the bottom panel
of Fig. <ref>, the time it takes for a fermion to travel
to the current measuring position, t_b, also marks shoulders
in the current (marked by the arrows), a much weaker effect than in
the partitioned case. As a matter of fact, in the partition-free case
weak features (inflection points or peaks) in the current can be observed
at both times t_b^L and t_b^R
as is shown in the bottom panel of Fig. <ref>.

After this initial build-up, in both approaches, the current enters
a sample-specific damped oscillatory phase which stabilizes towards
an approximately time-independent value. This stabilization marks
the beginning of the quasi-steady-state regime. As indicated in Fig. <ref>
and shown by the collapse of the curves in Fig. <ref>,
for different sample sizes, the quasisteady state exists for times
roughly greater than a stabilization time t_stab=2L_s/v_F.
Physically, this time can be interpreted as the one needed for a fermion
near the initial reference Fermi level to make a round trip inside
the central sample, thus probing the existing disorder landscape.
The fact that both t_b and t_stab are ballistic
times (i.e., ∝ L_meas and ∝ L_s,
respectively) is consistent with the fact that we are always working
in the ballistic regime of the mesoscopic central sample. Interestingly,
this time-scale is nearly independent of the particular disorder configuration
and the applied bias Δ V.











 §.§ Landauer Quasisteady state transport, finite-size effects and recurrence
times


In both approaches, if the leads are large enough, after the initial
build-up and stabilization of the current, a quasisteady state is
reached for t>t_stab, during which the current is approximately
time independent. As the size of the leads increases, the value of
this quasi-steady-state current tends to the sample-specific Landauer
value, independently of the initial preparation of the system (partitioned
or partition-free). Hence, the present results are numerical
checks to an extension of the memory-loss theorem of Stefanucci et
al. <cit.> for the case of finite
leads. The memory-loss theorem <cit.>
states that, provided the leads have a continuum spectrum, a steady-state
value of the current is achieved in the t→∞ limit,
and that this value is independent of the initial state of the system.

As can be seen in Figs. <ref> and
<ref>, in our case, and due to the finite nature
of the leads, which makes their spectrum discrete, a quasi-steady-state
only exists in a finite window of time: t_stab<t<t_r.
For t>t_r, we observe a drop and inversion of the current.
Similar behavior has also been observed previously <cit.>.
Pal et al <cit.> pointed out that the
recurrence time, t_r, is inversely proportional to the
level spacing of the leads' spectra, which measures how close the
finite leads are to a true continuous spectrum. Our results allow
for an alternative interpretation. As demonstrated in Figs. <ref>
— where we show the current for fixed Fermi energy and different
sizes of the leads — and in Fig. <ref> —
where we show the current for fixed L_l but different Fermi energies
— the recurrence time is roughly given by t_r=2L_l/v_F,
where v_F is the Fermi velocity. Notice that 2L_l/v_F
is just the time a fermion close to the Fermi level takes to perform
a round trip inside of a lead, in agreement with what was previously
reported in Ref. <cit.>. Furthermore, one
also sees that the recursion time is roughly independent of the disorder
on the sample, which is consistent with its previous physical interpretation.









One expects that, for sufficiently large leads, the value of the current
in the quasi-steady-state should approach the Landauer value, possibly
with some small oscillations superimposed due to finite size effects
that vanish with increasing lead size. Indeed, this is what we observe
in Fig. <ref>. However, the way
in which the oscillations vanish depends crucially on the initial
condition of the system. In the partition-free approach,
we observe that the current in the quasi-steady-state regime does
not strictly approaches the Landauer value. Instead, there is a small
persistent oscillatory component with a constant amplitude in time,
superposed on its steady-state value, and which persists up to the
recursion time t_r. Although constant throughout the
quasi-steady-state regime, the amplitude of these oscillations is
seen to decrease as L_l→∞, and the value of the quasi-steady-state
current approaches the Landauer value in a nearly uniform way for
t_stab<t<t_r. We also observed, that the period
of these oscillations is roughly inversely proportional to the applied
bias, i.e., T_osc∝Δ V^-1, but does not depend
on either L_l or L_s. Such dependence of the oscillation
period with the applied bias had also been reported by Kurth et
al <cit.>, although their focus is on systems with
infinite leads. Nevertheless, this behavior gives a hint on the physical
origin of these oscillations. If we consider an occupied eigenstate
|Ψ(0)⟩ of ℋ^PF(t=0), having
energy ε_0, this can always be written as a linear
combination of eigenstates |Ψ̃_n⟩ of ℋ(t>0).
This way, its time-evolved ket, |Ψ(t)⟩, is simply





    |Ψ(t)⟩=∑_ne^-iħ^-1ε_nt⟨Ψ̃_n|Ψ(0)⟩|Ψ̃_n⟩,

where ⟨Ψ̃_n|Ψ(0)⟩ are the wave-function
overlaps. Since the leads are much larger than the central sample
and the effect of applying a bias is to globally shift the energy
of all Wannier states in the leads to ±Δ V/2, it is reasonable
to say that the two dominant overlaps in Eq. (<ref>)
will be with the states |Ψ̃_±⟩ with energies
ε_0±Δ V/2. This follows from the fact that these
two are precisely the ones whose real-space wave-functions in either
the right or left lead are stationary waves with the same wavelength
as the original one. Assuming this argument to be true, we may neglect
all the other overlaps in Eq. (<ref>) and write





    |Ψ(t)⟩   ≃ e^-iħ^-1(ε_0+Δ V/2)t⟨Ψ̃_+|Ψ(0)⟩|Ψ̃_+⟩
        +e^-iħ^-1(ε_0-Δ V/2)t⟨Ψ̃_-|Ψ(0)⟩|Ψ̃_-⟩,

which behaves as a two level system, with a Bohr frequency ω=Δ V/ħ.
This argument justifies the presence of the time scale associated
to the period of the observed finite-size oscillations. Furthermore,
as Δ V is reduced, we expect that the overlaps ⟨Ψ̃_±|Ψ(0)⟩
will increase, thus, we also expect that the amplitude of these oscillations
will increase when Δ V decreases, a point to which we will
return.

For the partitioned approach, a rather different behavior is
observed. In this setup, as can be seen in the insets of Figs. <ref> (c)
and (d), the amplitude of the oscillations in the quasi-steady-state
decays as time increases, provided t_stab<t<t_r.
Furthermore, the amplitude of the oscillations is nearly independent
of the leads' size at any fixed observation time (provided t<t_r
for the compared lead sizes). As L_l increases, t_r=2L_l/v_F
also increases and therefore, the oscillations in the quasi-steady-state
will decay for a longer time, thus tending towards the Landauer value
as time tends to t_r, t→ t_r→∞.
Note that the physical argument given above in Eqs. (<ref>)
and (<ref>), justifies why one does not observe
persistent oscillations when the system begins in a partitioned setup.

In a true steady state, the value of the current is not only time
independent but must also be position independent, as no charge accumulation
can occur. Hence, we also investigated whether or not this emergent
quasisteady state current in finite chains is homogeneous over the
sample. Indeed, we found out that in the quasisteady state the current
is approximately homogeneous in space, independently of the initial
preparation of the system, for large enough leads and provided we
are far away from the chain's open extremities. This
observation is exemplified in Fig. <ref>, where
we show the time-dependent current for a disordered central sample,
measured at three different bonds: center, left, and right boundaries
of a randomly picked disordered sample. As can be seen, after the
disappearance of the initial transients, the same quasi-steady-state
current is reached at the three positions, apart from the finite-size
oscillations which are out of phase.









Interestingly, the establishment of a well-defined quasisteady state,
for a large but fixed leads size, might not occur for very small biases,
where we would expect linear response theory to hold, depending on
the initial conditions. This is illustrated in Fig. <ref>.
There, we can see that for the partition-free setup, no clear
quasi-steady-state is observed for very small biases. This occurs
because, for a fixed lead size and as previously explained, the period
and amplitude (relative to the infinite leads' Landauer value) of
the finite-size oscillations increases with the reduction of the applied
bias. Therefore, for small enough bias, the period of the oscillations
might become larger than the recursion time, and no quasi-steady-state
is observed. A well developed quasi-steady current only emerges provided
T_osc≪ t_r. In the partitioned setup,
the situation is a bit different and for large times: the current
always tends to the Landauer value with the amplitude of the finite-size
oscillations decreasing over time. These observations seem to be in
agreement with the interpretation of Bushong et al. <cit.>,
where it is put forward that the observation of a quasi-steady-state
requires the change in the initial spread of the electrons' momenta.
In this reference, this occurs either due to a geometrical constriction
at the lead-sample contact, or due to an initial applied energy barrier.
In our case, it seems that the applied bias is the mechanism by which
electrons change their initial momenta. As a general “rule-of-thumb”,
we can tell that, in order to observe a quasi-steady-state current
regime with minor finite-size effects, one must always consider biases
that are much larger than the level spacing of the whole system's
spectrum.






 §.§ Sample-Specific I-V Curves at Large Biases


We finally point out, that the coincidence between the Landauer value
for the current and the value of the current in the quasisteady state
occurs for any value of the bias potential, as long as it is smaller
than the bandwidth and provided a quasisteady state is established.

This is illustrated in Fig. <ref>, where we show values
for the time-dependent current in the quasisteady state regime as
a function of the applied bias, for two random disordered samples,
and compare the results with the value of the Landauer current. The
results clearly confirm that the quasisteady state current seen in
the quantum dynamics calculations with finite leads indeed corresponds
to the Landauer transport predicted for samples coupled to semi-infinite
leads. The agreement between the two approaches was seen to be perfect
for all the range of bias tested and well beyond linear response.









§ EMERGENCE OF LANDAUER TRANSPORT
WITHIN LINEAR RESPONSE IN THE PARTITION-FREE APPROACH


The numerical studies of the previous section show that a quasi-steady-state
transport regime, with an approximately uniform and time-independent
current, emerges across finite systems subjected to a potential bias
and coupled to finite but large leads. Moreover, the results also
show that for large enough leads, the value of this quasi-steady-state
current coincides with the Landauer result for the transport's steady
state with semi-infinite leads. In this section, we will try to shed
further light on these numerical results using a semi-analytical procedure.
In order to make as much analytical progress as possible, we shall
restrict ourselves to the partition-free case and small biases,
such that we can study the current using Kubo linear response theory
in the applied bias, Δ V.



 §.§ Time-dependent Kubo formula for a sudden connection


We will always consider the partition-free Hamiltonian at
t=0 as the unperturbed Hamiltonian for this case, i.e.,




    ℋ_0    =ℋ^PF(t=0)
        =∑_n=0^L-1ϵ_n^dc_n^†c_n-w∑_n=0^L-2(c_n+1^†c_n+c_n^†c_n+1)

and treat the applied potential bias as the current-driving perturbation,





    𝒱(t)=-eΘ(t)∑_n=0^L-1v_n^ec_n^†c_n.

with the electrostatic potential profile, v_n^e, given
by Eq. (<ref>).

In order to derive a time-dependent Kubo formula for the current,
we will start by writing the equation of motion for the reduced density
matrix, Eq. (<ref>), in the eigenbasis of the unperturbed
Hamiltonian. Thus, we obtain

    d/dtρ_αβ(t)=-i/ħ(ε_α-ε_β)ρ_αβ(t)-i/ħ[𝒱(t),ρ(t)]_αβ,

where O_αβ(t)=⟨Ψ_α|O|Ψ_β⟩
and |Ψ_α⟩ is an eigenstate of ℋ_0
with energy ε_α. Within linear response theory,
we write the reduced density matrix as

    ρ_αβ(t)=δ_αβf(ε_α)+δρ_αβ(t),

where ρ_αβ(0)=δ_αβf(ϵ_α)
is the initial equilibrium reduced density matrix and δρ_αβ(t)
is a small correction, which in linear response is assumed to be ∝𝒱(t).
Disregarding any contributions of 𝒪(𝒱^2)
in the equation of motion, we obtain





    d/dtδρ_αβ(t)=-i/ħ(ε_α-ε_β)δρ_αβ(t)
    
    -ie/ħΘ(t)Γ_αβ(f(ε_α)-f(ε_β)).

where Γ_αβ are the matrix elements of the applied
potential bias,




    Γ_αβ=∑_nψ_α^*(n)ψ_β(n)v_n^e,

and ψ_α(n) is the amplitude of the eigenstate |Ψ_α⟩
on site n, i.e. ψ_α(n)=⟨ n|Ψ_α.⟩.
Now, using the fact that δρ_αβ(t<0)=0, it is
possible to integrate Eq. (<ref>), obtaining

    δρ_αβ(t)=-eΓ_αβΔ f_αβ/Δε_αβ(1-e^-i/ħΔε_αβt),

where Δ f_αβ=f(ε_α)-f(ε_β)
and Δε_αβ=ε_α-ε_β.
The expected value of the current that flows from site n to n+1,
is thus given by

    I^n(t)=ie^2w/ħ∑_α,βΠ_αβ^nΓ_βαΔ f_αβ/Δε_αβ(1-e^-i/ħΔε_αβt),

where we introduced 

    Π_αβ^n=ψ_α^*(n+1)ψ_β(n)-ψ_α^*(n)ψ_β(n+1),

which are the matrix elements of the local current operator between
sites n and n+1, up to a dimension-full multiplicative factor.

By further noticing that the amplitudes ψ_α(n) may be
chosen as all real and Π^s is an anti-symmetric matrix, one
can rewrite Eq. (<ref>) in the following way:





    I^n(t)=2e^2w/ħ                              ∑_α(ε_α≤ε_F)∑_β(ε_β>ε_F)Π_αβ^nΓ_αβsin(Δε_αβt/ħ)/Δε_αβ.

which is our final time-dependent Kubo formula for the current.

Obviously, one cannot give a general rule for establishing the validity
regime of Eq. (<ref>), since that will depend crucially
on the properties of the central disordered sample. However, for each
sample, there is always a value of Δ V sufficiently small,
such that a linear response theory for the current is valid. We depict
such an example in the upper panel of Fig. <ref>,
where the current traversing the central bond of a disordered sample,
as obtained from Eq. (<ref>), is compared with
the one obtained from the fully nonlinear quantum dynamics of sec. <ref>
in the partition-free approach. As a further short comment
on the plots of Fig. <ref>, it is interesting
to note that, for the parameters used, it seems that no quasi-steady-state
plateau emerges from the quantum dynamics close to the linear response
regime. As referred before, this is simply a consequence of a greater
relevance of the finite-size oscillations which, now, have a period
larger than the recurrence time and a much larger relative amplitude.



 §.§ Representation of the eigenstates in terms of the sample's quantum
reflection/transmission coefficients


In order to make an effective use of Eq. (<ref>)
and make analytic progress we must be able to find a semi-analytical
expression for the matrix elements Π_αβ^n and Γ_αβ,
which, in principle, requires the knowledge of the eigenfunctions
in the whole chain. These wave functions usually present a very complicated
structure inside the disordered central sample, but for large enough
leads, we actually only need to know their form in the leads. On the
one hand, the Π_αβ^n matrix elements only require
the knowledge of local amplitudes in the two adjacent sites across
which the current is being measured. Hence, we can simply choose to
measure it outside the sample. On the other hand, we expect the current
to be dominated by states that are not localized in the disordered
sample, but instead are delocalized in the leads. Hence, we only need
to calculate the Γ_αβ matrix elements between delocalized
states. For such states, and provided the leads are much larger than
the disordered sample region, we can approximate


    Γ_αβ=∑_nψ_α^*(n)ψ_β(n)v_n^e≃    ∑_n∈Leadsψ_α^*(n)ψ_β(n)v_n^e.

This approximation, allows us to evaluate the current I^n(t)
in the leads, without knowing the shape of the eigenwave functions
inside the central sample.

Next, we notice that the form of the scattering eigenstates in the
leads can be expressed in terms of the complex reflection and transmission
coefficients of the central sample. For perfect leads, the wave functions
of the eigenstates will have the form of a coherent superposition
of left and right propagating plane-waves. With a change of notation
from the previous section, we will relabel sites of the left lead
with indices n=-L_l,...,-1 and the ones of the right lead with
n=1,...,L_l. Using this notation, the form of the eigenstate
wave function |Ψ_k⟩ in the leads have the
form





    ψ_k(n)=⟨ n|Ψ_k.⟩
    
    =Ψ_+^Le^ik(n-1)+Ψ_-^Le^-ik(n-1),              -L_l≤ n≤-1
    Ψ_+^Re^ik(n-L_s)+Ψ_-^Re^-ik(n-L_s),     1≤ n≤ L_l,

being labeled by a crystal momentum k, and with Ψ_+/-^L(R)
being the amplitude of a right/left propagating state in the left
(right) lead. Notice that the time-independent Schrdinger equation
inside the leads, still allows us to relate the crystal momentum k
to the energy of the state as E=-2tcos(k), i.e. the same as for
an infinite periodic chain. As usual in one-dimensional scattering
problems, the amplitudes of propagating states on the left and right
leads can be related by a transfer matrix, ℳ(k):

    ([ Ψ_+^R; Ψ_-^R ])=ℳ(k)·([ Ψ_+^L; Ψ_-^L ]),

In the presence of time-reversal symmetry, the transfer matrix has
the general form

    ℳ(k)=([            1/t(k)e^iϕ(k) -r(k)/t(k)e^-iθ(k)+iϕ(k);  -r(k)/t(k)e^iθ(k)-iϕ(k)           1/t(k)e^-iϕ(k) ]),

where t(k)/r(k) and ϕ(k)/θ(k)
are the moduli and phases of the transmission and reflection coefficients,
respectively. Moreover, for any sample one has ℳ=1,
which implies the conservation of current, i.e., t^2+ r^2=1.
These coefficients are physical characteristics of the central sample
only and, thus, may be rightfully calculated by assuming the leads
as semi-infinite. The determination of the reflection and transmission
coefficients of a specific sample, in general, can only be done numerically,
using the method detailed in Appendix <ref>.
The great advantage of this method is that, once this calculation
is done, the wave functions in the leads can be expressed in terms
of only a few parameters. Additionally, to obtain the eigenstates,
we must further impose open boundary conditions at the ends of the
leads, i.e.,




    ψ_k(-L_l-1)=ψ_k(L_l+1)=0.







Combining Eqs. (<ref>)-(<ref>)
one arrives at the following general expression for the wave functions:





    ψ_k(n)=1/√(N_k)|t(k)|sin[k(n+L_l+1)]     n<0
    
    f_2(k)sin[k(n-L_l-1)]     n>0
    ,

where N_k is a normalization factor, which can be determined
in the limit of large leads by approximating, in the same spirit of
Eq. (<ref>),

    ∑_n|ψ_k(n)|^2≃∑_n∈leads|ψ_k(n)|^2.


This finally leads to




    N_k≃ L_lf_1(k).

The functions f_1(k) and f_2(k) are defined as







    f_1(k)     =1+r(k)cos[2k(L_l+1)+θ(k)]
    
    f_2(k)     =cos[2k(L_l+1)+ϕ(k)]
        +r(k)cos[θ(k)-ϕ(k)],

 and where k is constrained to verify the following
quantization condition:



    sin[2k(L_l+1)+ϕ(k)]=r(k)sin[θ(k)-ϕ(k)].



Notice that the solution of this last condition, together with the
relation k=arccos(-E/(2w)), allows us to
determine the eigenenergies corresponding to delocalized states. In
the lower panels of Fig. <ref> (b),
we exemplify the validity of this statement by comparing the wave
functions obtained from the numerical diagonalization of ℋ_0,
with L_l=8192 sites, to the semi-analytical expressions of Eq. (<ref>).
The wave numbers obtained from the numerical diagonalization, i.e.
k=arccos(-E/2), are also seen to coincide perfectly
with the roots of Eq. (<ref>) [see the
upper panels of Fig. <ref> (b)].

Note also that Eqs. (<ref>) and (<ref>)
reduce to the usual result for the eigenstates of a finite open chain,
when r(k)=0 and ϕ(k)=k(L_l-1)
is the phase accumulated by a plane-wave crossing the internal bonds
of an ordered sample, i.e.,


    ψ_k(n)=1/√(L_l)sin[k(n+L_l+1)],      n<0
    (-1)^psin[k(n-L_l-1)],n>0

with k=π p/(L+1) and p=1,...,L. Moreover, these
states are non-degenerate and also alternately symmetrical and antissymmetrical
under parity (n→-n), which just reflects that same symmetry of
the clean Hamiltonian.

With the knowledge of the eigenstates wave functions of the leads,
Eq. (<ref>), we can write the matrix elements
Γ_k,q and Π_k,q^n. With the approximation of Eq. (<ref>),
we can evaluate Γ_k,q analytically, obtaining

    Γ_k,q≃Δ Vt(k)t(q)-f_2(k)f_2(q)/8L_l√(f_1(k)f_1(q))×
    ×{sin[(q-k)(L_l+1/2)]/sin(q-k/2)-sin[(q+k)(L_l+1/2)]/sin(q+k/2)} .



As for the matrix elements Π_k,q^n, from the definition Eq. (<ref>),
for bonds in the left lead (n<-1), and after some simple manipulations,
we obtain

    Π_k,q^n<-1=|t(k)||t(q)|/L_l√(f_1(k)f_1(q))×
    ×{sin(k-q/2)sin[(k+q)(n+L_l+3/2)].
    .-sin(k+q/2)sin[(k-q)(n+L_l+3/2)]} ,

while for the current in bonds of the right lead (n>1), we obtain
a similar result after replacing |t(k)||t(q)|→ f_2(k)f_2(q)
and L_l→-L_l-2 in Eq. (<ref>).



 §.§ Continuum regime of the Kubo formula


When analyzing the time-dependent Kubo formula of Eq. (<ref>),
we must take into account that there are actually two distinct time
scales: (1) the observation time, t, and (2)
the scale associated with the spacing between the discrete energy
levels of the finite chain. The latter is proportional to the length
of the leads and, as discussed in the Sec. <ref>,
is associated to the recurrence time t_r∼2L_l/v_F.

As expected and confirmed in Sec. <ref>, the
quasi-steady-state regime which approximates the Landauer transport
regime of semi-infinite leads, emerges when we take T,L_l→∞
(with T=tw/ħ being the time in dimensionless units), but while
keeping T≪ L_l. In such case, all the transients have died
out, but the system is still far away from getting into the regime
where current inversions occur. Furthermore, Eq. (<ref>)
includes a factor of sin(Δε_αβt/ħ)/Δε_αβ,
which is an emergent δ-function in the limit t→+∞,
with a broadening of ħ t^-1 in energy. This factor actually
acts as a spectral filter which kills-off the contributions coming
from pairs of eigenstates having an energy separation larger than
ħ t^-1. Hence, we will show in this section how the approximately
time-independent quasi-steady-state current emerges, when we are in
the limit T,L_l→∞, with T≪ L_l, such that
there are many eigenvalues inside the interval [ε_F-ħ t^-1,ε_F+ħ t^-1].
We will refer to this limit as the continuum regime.






  §.§.§ Approximate form of Γ_k,q and Π_k,q^n matrices
in the continuum regime


We start by noting that, in the continuum limit, since only states
close to the Fermi energy contribute, it suffices to obtain the matrix
elements Γ_k,q and Π_k,q^n between states where
k-q is small and k,q≃ k_F. In the limit of k-q→0,
the first term of Eq. (<ref>) dominates over
the second. Therefore, we can approximate it as




    Γ_k,q≃Δ Vt(k)t(q)-f_2(k)f_2(q)/8L_l√(f_1(k)f_1(q))sin[(k-q)L_l]/sin(k-q/2),

where, in the of limit L_l→∞ , we approximated
sin[(q-k)(L_l+1/2)]≃sin[(q-k)L_l].
Doing the same for Π_k,q^n, we obtain




    Π_k,q^n<-1≃-t(k)t(q)/L_l√(f_1(k)f_1(q))sin(k_F)sin[(k-q)L_l],

where we assumed that |n|≪ L_l, when approximating
sin[(k-q)(n+L_l+3/2)]≃sin[(k-q)L_l].
This justifies why in the quasi-steady-state regime, the current is
approximately uniform, if we are away from the chain's extremities.
For the current on the right lead, we obtain a similar result, namely,




    Π_k,q^n>1≃f_2(k)f_2(q)/L_l√(f_1(k)f_1(q))sin(k_F)sin[(k-q)L_l].


Now, we note that for a chain without any disorder, the matrix elements
of Γ_k,q will only be non-zero if the states labeled to
k and q have opposite parities. This selection rule stems from
the fact that the fully ordered chain is symmetric under inversion
and therefore its eigenstates will have a well-defined parity. Since
the applied potential v_n^e is an odd perturbation,
it only couples states of opposite parities. In the presence of a
general disorder in the central sample, we no longer have inversion
symmetry. Nevertheless, one may still expect that in the limit L_l≫ L_s,
the breaking of the symmetry is small and an approximate selection
rule should emerge. Indeed, this is the case. In order to obtain this
approximate selection rule for a sample with disorder, we notice that
although we can no longer classify the states as even and odd, given
the quantization condition Eq. (<ref>),
which involves sin[2k(L_l+1)+ϕ(k)],
we can classify the states as + and - according to the sign
of cos[2k(L_l+1)+ϕ(k)]:

    cos[2k^±(L_l+1)+ϕ(k^±)]=
    
    =±√(1-r(k^±)^2sin^2[θ(k^±)-ϕ(k^±)]).





For an ordered or symmetrically disordered sample, this
reduces to a labeling of eigenstates as even or odd, respectively,
under a parity transformation, n→-n. With such a classification,
it can be shown (see Appendix <ref>) that in
the limits of k-q→0 and L_l→∞, one
obtains the following effective selection rule:

    lim_L_l→∞|sin[(k^σ-q^σ^')L_l]|=(1-δ_σ,σ^')×
    ×√(1-r(k_F)^2sin^2(θ(k_F)-ϕ(k_F))).




with σ,σ'=± and which immediately implies
that Γ_k,q≃0, if k,q are in the same class as L_l→∞.
This is an approximate selection rule, analogous to the one which
exists in the clean case, but which only emerges when L_l→∞.
In Fig. <ref>, we represent the values of |sin[(k-q)L_l]|
as a function of (k-q)L_l, for allowed values of k
and q. We can clearly see that for some data points |sin[(k-q)L_l]|→0
as L_l increases, while other data points tend to a finite value,
which is given by the sample-specific value, √(1-r(k_F)^2sin^2(θ(k_F)-ϕ(k_F))). [As visible in Fig. <ref>, there is small deviation
of the data points obtained for the ordered central sample from the
theoretical value sin(L_l(q^σ-k^σ'))=1-δ_σ,σ'.
These results seem incompatible with an exact parity selection rule
for Γ_k,q, for this case, however, they are not. This artifact
is due to the fact that the data shown was calculated from the numerical
diagonalization of ℋ_0 and then k=arccos(-E/2)
was used to obtain the respective wave numbers. This procedure takes
into account the finite dimension of the central sample and hence
the allowed wave numbers are of the form k^σ=π n/(2L_l+L_s+1),
with σ=(-1)^n. Since the expression of Γ_k,q
which is proportional to sin(L_l(q-k))
is only valid under an approximation which ignores L_s [i.e.
Eq. (<ref>)], the parity selection rule
appears to be approximate as well. At any rate, it may be proven,
by symmetry, that this rule is actually true for any value of L_l,
if one considers the full expression for Γ_k,q.]

We point out that in the case of symmetric disorder profile, one can
derive from the properties of the transfer matrix that ϕ(k)-θ(k)=±π/2.
In this case, one immediately sees that the scattering wave functions
of Eq. (<ref>) reduce to the same form as
in Eq. (<ref>), with the parity determined
by the class to which it belongs. In such a case, the t(k)
factor of the Γ and Π matrices comes only from this effect,
since the functions f_1(k) and f_2(k)
are exactly the same as in the non-disordered case. Just the allowed
k's are different.

Having established this effective selection rule, we can expand the
prefactors of Eqs. (<ref>)—(<ref>)
around k_F. Taking into account that the only significant
contributions come from pairs of states belonging to different classes,
we may use Eqs. (<ref>) and (<ref>)
to write, for k,q≃ k_F,






    f_1(k)f_1(q)≃
    ≃(1-r(k_F)^2sin^2[θ(k_F)-ϕ(k_F)])|t(k_F)|^2,


    f_2(k)f_2(q)≃-|t(k_F)|^2.

Using these approximations in Eq. (<ref>),
we obtain

    Γ_k,q   ≃Δ V|t(k_F)|/4L_l1/sin(k-q/2)
       ≃Δ V|t(k_F)|/2L_l1/k-q,

and from Eqs. (<ref>) and (<ref>)
we obtain

    Π_k,q^n<-1≃Π_k,q^n>1≃-|t(k_F)|/L_lsin(k_F).

In the following, we will use Eqs. (<ref>)
and (<ref>) to obtain the Landauer current from
the Kubo formula of Eq. (<ref>)





  §.§.§ Continuum limit expression for the stationary current: Emergence
of Landauer transport


Using Eqs. (<ref>) and (<ref>),
we can write the time-dependent Kubo formula Eq. (<ref>)
as

    I^n(t)=e/2^2^'∑_k,q𝒫_k,qsin(Δε_k,qt/ħ)/Δε_k,q,

where the primed sum in Eq. (<ref>) means that
only pairs of states (k,q) of opposite classes are included
in the sum, due to the emergent selection rule. We also introduced
the quantity 𝒫_k,q, which is defined as

    𝒫_k,q    =2w/ħΠ_k,q^nΓ_k,qΔ f_k,q
       ≃-1/2L_l^2|t(k_F)|^2v_F^2ħ/Δε_k,qΔ V,

where we approximated v_Fħ(k-q)≃Δε_k,q,
with v_F=2wsin(k_F)/ħ. To make
Eqs. (<ref>) and (<ref>)
more clear, we remark that this definition of the current is no longer
dependent on the condition k<k_F<q, and 𝒫_k,q
is actually symmetrical upon exchange of the indices. The above equation
also shows that, for k,q≃ k_F, the latter approximately
only a function of the difference in eigenenergies. This result was
checked numerically as seen in Fig. <ref>, where
we can see that for a wide variety of disordered samples and different
values of the Fermi energy, all the values of 𝒫_k,q
[calculated directly from the wave functions in Eq. (<ref>)]
fall into the curve Eq. (<ref>).

The time-dependent current in the continuum regime can thus be written
as





    I^n(t)=e^2/ħ|t(k_F)|^2(v_Fħ)^2Δ V
    ×∫_-∞^+∞d(Δε)sin(Δε t/ħ)/ΔεΔεϱ(Δε),

where we introduced the joint density of contributing states
(JDoCS), ϱ, as





    ϱ(ϵ_F,Δε)=1/4L_l^2^'∑_k,qδ(Δε-Δε_k,q).


The restricted summation in Eq. (<ref>) already takes
into account the emergent selection rule of Eq. (<ref>).
In Appendix <ref>, we show that this quantity,
in the limit L_l→∞, can be written in terms of the density
of states of each class in a fully clean system and its expression
for small enough Δε is simply





    lim_L_l→∞[ϱ(ε_F,Δε)]=Δε/2π^2(4w^2-ε_F^2)+𝒪[Δε^2].

Hence, when Eq. (<ref>) is plugged into Eq. (<ref>),
we get


    I^n(t)=e^2/2π^2ħ|t(k_F)|^2(v_Fħ)^2Δ V/4w^2-ε_F^2
    ×∫_-∞^∞d(Δε)sin(Δε t/ħ)/Δε.

Finally, Eq. (<ref>) together with the facts
that





    lim_T→∞[sin[xT]/x]=πδ(x),

and v_Fħ=√(4w^2-ε_F^2), yields a steady-state
current

    I^n(t)=e^2/h|t(k_F)|^2Δ V,

which is precisely the linear Landauer steady-state current for a
two-terminal one-dimensional device.

Notice, that in the derivation of this result from the time-dependent
Kubo formula, it is essential that t,L_l→∞ with
wt/ħ≪ L_l, such that the ϱ(Δε)
can be evaluated in the limit of L_l→∞, while
the factor sin(Δε t/ħ)/Δε
is treated as as emergent δ-function. When wt/ħ L_l,
then there will be few pairs of states with Δε_k,q∈[ε_F-ħ t^-1,ε_F+ħ t^-1],
and we can no longer treat sin(Δε t/ħ)/Δε
as a δ-function. When this happens, we start observing recurrences
in the current as reported in Sec. (<ref>).






§ CONCLUSIONS


In this work, we investigated how a quasi-steady-state particle transport
regime emerges across disordered samples coupled to large, but finite
leads which are subjected to a potential bias. In order to do so,
we have studied time-dependent transport, both numerically and semi-analytically,
in a non-interacting and one-dimensional tight-binding chain, with
open boundary conditions, where the central region is an extended
disordered sample, and the rest of the chain acts as a pair of finite,
but otherwise perfect leads.

For large lead sizes, and sufficiently large bias, a quasi-steady-state
regime emerges at intermediate times, after the transient behavior
has died out and before inversions in the current are observed. The
current in the quasi-steady-state is approximately constant in time
and homogeneous in space (if measured at points far away from the
chain's extremities). Furthermore, the value of the current in the
quasi-steady-state coincides with the one predicted by the Landauer
formula for semi-infinite leads, independently of the initial condition
of the system (partitioned or partition-free). These
results amount to an exemplification and extension to finite systems
of the results of Stefanucci et al <cit.>
on the establishment of a steady-state regime of transport in samples
which are attached to infinite leads.

We have found that the quasisteady state is established, for both
initial conditions, after a stabilization time t_stab≈2L_s/v_F.
Physically, this can be interpreted as the time taken by a Fermi-level
state to probe the disordered landscape inside the central sample.
The quasi-steady-state lasts until a recurrence time t_r≈2L_l/v_F,
where current inversions start happening. Aside from being related
to the inverse spacing of the energy levels in the system <cit.>,
this recurrence time may also be interpreted as the time taken by
a Fermi-level electron to leave the sample and return to it, by traveling
back and forth inside a lead. This conclusion was seen to be independent
of the central sample's features, as long as the leads are much larger
than it and transport is ballistic across the disordered sample.

During the quasi-steady-state, persistent finite-size effects are
observed in the partition-free approach as superposed oscillations,
with a period that is inversely proportional to the bias Δ V
and an amplitude which scales to zero as L_l→∞
but becomes more relevant (relative to I_Landauer) for
very small values of Δ V. This effect prevents the onset of
a quasi-steady-state regime for systems prepared in the partition-free
setup, if the leads are too small. In the partitioned case,
the amplitude of the oscillations superposed on the quasisteady-state
plateaus is not influenced by the size of leads, but instead is damped
as the observation time increases (while keeping t<t_r).
Similarly to the partition-free case, the amplitude of the
fluctuations increases for smaller biases. These observations seem
to indicate that the observation of a clear quasisteady state requires
some kind of mechanism which scatters the electron's momenta <cit.>.
Here it is provided by the applied potential ramp in the sample, which
becomes a less effective mechanism as Δ V→0. In both cases,
these finite-size oscillations can be made arbitrarily small if L_l
is large enough.

In order to shine light on the numerical results, a time-dependent
Kubo formula for the current in the partition-free approach,
which is suitable for semi-analytical treatment, was developed for
describing the local time-dependent current due to a small applied
bias. From this formula, it was possible to see that an approximately
time-independent and spatially uniform current emerges in the limit
of large system's size and observation times, L_l,t→∞,
provided t≪ v_F^-1L_l (in agreement with the recurrence
times observed numerically). These conditions are necessary to treat
the leads as being effectively infinite, in what respects DC transport.
After expressing the eigenfunctions of the disordered central sample
in terms of complex reflection and transmission coefficients, all
the matrix elements appearing in the Kubo formula were evaluated semi-analytically.
The quasi-steady-state current thus obtained was shown to reproduce
the linearized Landauer formula for the current in a two-terminal
device.

We hope that these theoretical predictions of the time scales over
which the quasisteady state occurs and the nature of the finite-size
oscillations can be experimentally tested and guide future research
on mesoscopic transport in fermionic ultra-cold atomic gases in optical
lattices.






§ ACKNOWLEDGMENTS


J.M.V.P.L. and J.P.S.P. acknowledge financing of Fundao da Cincia
e Tecnologia, of COMPETE 2020 program in FEDER component (European
Union), through projects POCI-01-0145-FEDER-028887 and UID/FIS/04650/2013.
J.P.S.P. is supported by the MAP-fis PhD grant PD/BD/142774/2018 of
Fundao da Cincia e Tecnologia. B.A. acknowledges financial support
from Fundao para a Cincia e a Tecnologia, Portugal, through Project
N CEECIND/02936/2017. Additionally, J.P.S.P. also acknowledges the
hospitality of the University of Central Florida, where part of this
work was done, as well as Dr. Eduardo Mucciolo and Dr. Caio H. Lewenkopf
for the careful reading of the manuscript and suggestions. We also
thank S. M. Joo, Dr. Nuno M. R. Peres, Dr. Aires Ferreira and Dr.
Joo M. B. Lopes dos Santos for enlightening discussions about this
work. J.P.S.P. also acknowledges the hospitality of Katherine Vasquez
during his stay in Orlando, Florida. Finally, the authors acknowledge
the comments of the anonymous referees, which proved very helpful
in the improvement of the original manuscript.





§ REVIEW OF THE RECURSIVE CHEBYSHEV
METHOD FOR QUANTUM TIME-EVOLUTION


In this appendix, we wish to describe shortly the algorithm used to
time-evolve an arbitrary single-particle state with the full Hamiltonian.
As referred in the main text, the Hamiltonian generating the time-evolution
for positive times, ℋ(t>0), is time-independent
and, consequently, the time-evolution operator 𝒰_t
reads





    𝒰_t=e^-iℋ(t>0)t/ħ.


The method used to calculate 𝒰_t for our systems is
based on its exact expansion as a series of Chebyshev polynomials
in ℋ(t>0), due to Tal-Ezer et al.
<cit.>. Namely, one has





    𝒰_t=∑_n=0^∞2/1+δ_n,0(-i)^nJ_n(λ t)T_n(ℋ̃),

where ℋ̃=(1/λ w)ℋ
is a dimensionless Hamiltonian, rescaled by a real parameter λ
which guarantees that its spectrum is contained inside the interval
]-1,1[, T_n is the n^th-order Chebyshev
polynomial of the first-kind, J_n(y) is a Bessel function of
the first kind and t is a time measured in units of ħ/w.
The key to the method is to avoid the numerical diagonalization of
ℋ(t>0), and instead use the recursion relation
for the Chebyshev polynomials,





    T_n+1(x)=xT_n(x)-T_n-1(x),

in order to evaluate all the needed T_n(ℋ̃),
recursively. For a generic review on the application of Chebyshev
spectral method to physical problems see Ref. <cit.>
and references within.



Furthermore, the Chebyshev series of Eq. (<ref>)
is known to converge rather quickly, meaning that a truncated summation
with M terms is usually enough to describe correctly 𝒰_t,
provided M>tλ. This convergence is illustrated in Fig. <ref>
and in all our calculations, we used M=8tλ.

Notice that, in order to evaluate the current, we only require to
time-evolve a given single-particle state |Ψ⟩. Therefore,
we do need the full matrix form of 𝒰_t, but instead
how it acts on an arbitrary state |Ψ⟩. From the expansion
of Eq. (<ref>), we know that action to
be





    |Ψ^M(t)⟩=∑_n=0^M2/1+δ_n,0(-i)^nJ_n(λ t)|Ψ_n⟩,

where |Ψ_n⟩=T_n(ℋ̃)|Ψ⟩
and M is the truncation order of the Chebyshev expansion. [As this is needed in the main text, we remark that for backward time
evolutions, one may simply use the fact that J_n(-x)=(-1)^nJ_n(x).] Finally, the first two |Ψ_n⟩ can be directly calculated
by the simple forms of T_0(x) and T_1(x),
i.e.







    |Ψ_0⟩    =T_0(ℋ̃)|Ψ⟩=|Ψ⟩
    |Ψ_1⟩    =T_1(ℋ̃)|Ψ⟩=ℋ̃|Ψ⟩,

and then the remaining are efficiently calculated
by using the operator generalization of the Chebyshev recursion [Eq. (<ref>)],
i.e.





    |Ψ_n+1⟩=ℋ̃|Ψ_n⟩-|Ψ_n-1⟩.





§ REVIEW OF THE RECURSIVE
TRANSFER MATRIX METHOD


In this appendix, we explore a very simple algorithm which allows
us to calculate the transfer matrix ℳ(k) of
any given disordered sample, when it is connected to semi-infinite
leads. This method is the same used in the early papers of Andereck
et al <cit.> and Pichard <cit.>
and allows for the calculation of ℳ(k) with
an ∼𝒪(L_S) number of operations.

For these purposes, it is more useful to re-express the Hamiltonian
of the central sample in a first-quantization language, i.e.,




    ℋ_s    =∑_n=1^L_sε_n|n⟩⟨ n|
               -∑_n=1^L_s-1(|n⟩⟨ n+1|+|n+1⟩⟨ n|),

where |n⟩ are the Wannier states of the chain
and ε_n is an on-site energy (in units of the hopping
w). To model the connection between the finite sample to the semi-infinite
leads, one has also the following boundary hopping Hamiltonian:






    ℋ_s=-|0⟩⟨1|-|1⟩⟨0|-|L_s⟩⟨L_s+1|-|L_s+1⟩⟨L_s|.



The main purpose of this method is to find the scattering states associated
to a particular disorder realization. For that, one must fix the leads'
propagating states, |Ψ_±^L⟩ and |Ψ_±^R⟩,
as the left and right boundary conditions for the problem. This setup
is represented in Fig. <ref>, with the counter-propagating
plane waves in the leads being represented as arrows.








 §.§ Hamiltonian in Real-Space and Boundary Conditions


The first step towards the definition of the present method is expanding
a scattering state (with wavenumber k) in the basis of Wannier
wave functions, i.e.





    |Ψ_k⟩=∑_nψ_n|n⟩,

and finally rewriting the time-independent Schrdinger equation, ℋ|Ψ_k⟩=E_k|Ψ_k⟩,
in terms of the real-space amplitudes ψ_n,





    E_kψ_n=ε_nψ_n-ψ_n-1-ψ_n+1,

where, by definition, ε_n=0 outside of the sample.

As shown in Fig. <ref>, the boundary conditions
are to be set as the plane waves defined in Eq. (<ref>)
of the main text. Reminding, one has






    |Ψ_k^L⟩=∑_n=-L_l^-1[Ψ_+^Le^ik(n-1)|n⟩+Ψ_-^Le^-ik(n-1)|n⟩],








    |Ψ_k^R⟩=∑_n=1^L_l[Ψ_+^Re^ik(n-L_s)|n⟩+Ψ_-^Re^-ik(n-L_s)|n⟩].

These states immediately set the amplitudes on the “ghost” sites
of Fig. <ref> to the following values:





    ψ_-1    =Ψ_+^Le^-2ik+Ψ_-^Le^2ik,
    ψ_0    =Ψ_+^Le^-ik+Ψ_-^Le^ik,
    ψ_L_S+1    =Ψ_+^Re^ik+Ψ_-^Re^-ik,
    ψ_L_S+2    =Ψ_+^Re^2ik+Ψ_-^Re^-2ik.







 §.§ Review of the transfer matrix recursive method


Despite not having the look of a linear algebra problem, Eq. (<ref>)
may be turned into a matrix recursion equation, when supplemented
by the trivial condition





    ψ_n=ψ_n.

Hence, we have





    ([ ψ_n+1;   ψ_n ])=𝕋_n(k)([ ε_n-E_k      -1;       1       0 ])·([   ψ_n; ψ_n-1 ]).

If we now iterate Eq. (<ref>), we get the
following relation





    ([ ψ_L_s+2; ψ_L_s+1 ])=𝕋_L_s+1(k)·𝕋_L_s(k)·
    ⋯·𝕋_1(k)·𝕋_0(k)·([  ψ_0; ψ_-1 ]).


In the same way, we may write the boundary conditions of Eqs. (<ref>),
as the following matrix relations:





    [ ([  ψ_0; ψ_-1 ])=𝔹_L(k)([  e^-ik   e^ik; e^-2ik  e^2ik ])·([ Ψ_+^L; Ψ_-^L ]),;                                                                              ;                                                                               ]


and





    ([ ψ_L_s+2; ψ_L_s+1 ])=([  e^2ik e^-2ik;   e^ik  e^-ik ])·([ Ψ_+^R; Ψ_-^R ]),

which can be inverted as




    ([ Ψ_+^R; Ψ_-^R ])=𝔹_R(k)·([ ψ_L_s+2; ψ_L_s+1 ])


Using Eqs. (<ref>) and (<ref>) into
Eq. (<ref>), we get to the following final
result:





    ([ Ψ_+^R; Ψ_-^R ])=𝔹_R(k)·𝕋_L_s+1(k)·𝕋_L_s(k)·
    ⋯·𝕋_1(k)·𝕋_0(k)·𝔹_L(k)·([ Ψ_+^L; Ψ_-^L ]),

and, by definition, the transfer matrix of the whole sample is written
as:





    ℳ(k)=𝔹_R(k)·𝕋_L_s+1(k)·𝕋_L_s(k)·
    ⋯·𝕋_1(k)·𝕋_0(k)·𝔹_L(k).

This last equation was the one we implemented to calculate ℳ(k)
for any given disordered sample.





§ EMERGENCE OF SELECTION RULE


In this appendix, we prove the effective selection rule of Eq. (<ref>).
In order to do so, we will analyze the factor sin[(k-q)L_l],
when q,k belong to the same or different classes. More precisely,
will calculate its absolute value, which can be written as





    sin[(k-q)L_l]=√(1-cos[2L_l(q-k)]2)
    
    =1/√(2){ 1-cos[2(L_l+1)(q-k)+ϕ(q)-ϕ(k)]cos[ϕ(q)-ϕ(k)-2(q-k)].
    
    -.sin[2(L_l+1)(q-k)+ϕ(q)-ϕ(k)]sin[ϕ(q)-ϕ(k)-2(q-k)]} ^1/2,

where we summed and subtracted ϕ(q)-ϕ(k)
in the argument of the cosine and, then, decomposed it using the rule
for the cosine of a sum of angles. The main advantage of this form
is that the continuous function ϕ(k) depends solely
in the properties of the central sample and the effect of increasing
the leads is to populate more densely their domains with allowed values
of k. This, together with the fact that we are only interested
in what happens near k_F, allows us to expand it as Taylor
series on δ q=q-k_F and δ k=k_F-k:


    ϕ(q)-ϕ(k)=.d/dkϕ|_k_F(δ q+δ k)+⋯≃.d/dkϕ|_k_F(q-k),

and, consequently,





    sin[(q-k)L_l]≃   √(1-cos[2L_l(q-k)+ϕ(q)-ϕ(k)]2),

where the corrections are of order q-k and disappear in the limits
L_l→∞ and ħ t^-1→0. At this point, all we must
do is to decompose the cosine term in Eq. (<ref>) using
the usual rules for the sum of angles and then resort to the quantization
condition of Eq. (<ref>) to realize that





    cos[2(L_l+1)(q-k)+ϕ(q)-ϕ(k)]
    
    =∓√([1-r(q)^2sin^2(θ(q)-ϕ(q))][1-r(k)^2sin^2(θ(k)-ϕ(k))])
    
    +r(q)r(k)sin[θ(q)-ϕ(q)]sin[θ(k)-ϕ(k)],

where the +(-) sign stands for the case when q and
k are in the same class (different classes) of states.

Finally, one can evoke the same argument as before to Taylor expand
all the sample-specific functions appear in Eq. (<ref>)
[to be clear, r(x),θ(x),ϕ(x)]
around k_F, but noting that k<k_F<q by definition.
Up to corrections irrelevant correction in the same limits, this gives
rise to Eq. (<ref>) of the main text after
expanding the sin functions in powers of q-k.







§ CALCULATION OF THE JOINT DENSITY OF
CONTRIBUTING STATES


In this appendix, we will proceed to calculate the joint density of
contributing states (JDoCS), for both positive and negative Δε.
For positive energy differences, Δε>0, the JDoCS
is defined, from Eq. (<ref>), as





    ϱ(ϵ_F,Δε)=1/4L_l^2^'∑_k,q(ε_q<ε_F≤ε_k)δ(Δε-Δε_k,q),

which may be written in terms of the usual density of states for each
class, σ=±, i.e.,




    ρ^σ(ε)     =1/L_l∑_k^σδ(ε-ε_k^σ),

yielding the expression,




    ϱ(ε_F,Δε)=1/4∫_ε_F^2dϵ_2∫_-2^ε_Fdϵ_1lim_L_l→∞{ρ_L_l^+(ϵ_1)ρ_L_l^-(ϵ_2).
    .+ρ_L_l^-(ϵ_1)ρ_L_l^+(ϵ_2)}δ(Δε-ϵ_2+ϵ_1),

in the limit of semi-infinite leads.

To progress beyond Eq. (<ref>) in a general fashion, one
starts by recognizing that, since ρ^±(ε)
is an intensive quantity. So these must be dominated by the states
on the (clean) leads, as L_l→∞. Since we know that, for
a clean system, the states of different parities are alternated in
k-space, with a regular separation given by π/L_l, one concludes
that





    lim_L_l→∞ρ_L_l^±(ε)=ρ(ε)=1/π√(4w^2-ε^2)   ifε≤2w
    
    0    ifε>2w
    ,


where ρ(ε) is the full DoS of
a clean infinite chain.

In what follows, we will always assume that the expression of Eq. (<ref>)
may be used to calculate de JDoCS in the limit of very large L_l.
This intuition is confirmed by the plots of the DoS in Fig. <ref>,
which were obtained numerically, for a randomly selected disordered
sample, using the well-known kernel polynomial method with
a Jackson kernel and a fixed number of polynomials, M=4096, enough
to resolve the individual energy levels in the smaller case considered
(see Weie et al. <cit.> for more details
on the method). Consequently, one has the following expression for
the JDoCS:




    ϱ(ε_F,Δε)     =∫_ε_F^2dε_2Θ(Δε+ε_F-ε_2)/2π^2√((4w^2-ε_2^2)(4w^2-(ε_2+Δε)^2))
        =∫_ε_F^ε_F+Δε                            dε_2Θ(Δε+ε_F-ε_2)/2π^2√((4w^2-ε_2^2)(4w^2-(ε_2+Δε)^2))

where Θ(x) is the Heaviside function and Δε≥0.
The integral in Eq. (<ref>) can be done numerically
and the curves are shown in Fig. <ref> for different values
of the Fermi energy ϵ_F. Nevertheless, we are only
interested in the shape of ϱ(ε,Δε)
when Δε≈0. For that, we may expand Eq. (<ref>)
in powers of this quantity, yielding





    ϱ(ε_F,Δε>0)=Δε/2π^2(4w^2-ε_F^2)+𝒪[Δε^2].


Finally, we can generalize Eq. (<ref>) to Δε<0,
which is trivial since, by definition [Eq. (<ref>)],
we have ϱ(Δε)=ϱ(-Δε).
Hence, our final expression is simply,





    ϱ(ϵ_F,Δε)=Δε/2π^2(4w^2-ε_F^2)+𝒪[Δϵ^2],

which is the one we use in the main text [see Eq. (<ref>)].

















apsrev4-1


