Non-equilibrium quantum transport in presence of a defect: the non-interacting case

We study quantum transport after an inhomogeneous quantum quench in a free fermion lattice system in the presence of a localised defect. Using a new rigorous analytical approach for the calculation of large time and distance asymptotics of physical observables, we derive the exact profiles of particle density and current. Our analysis shows that the predictions of a semiclassical approach that has been extensively applied in similar problems match exactly with the correct asymptotics, except for possible finite distance corrections close to the defect. We generalise our formulas to an arbitrary non-interacting particle-conserving defect, expressing them in terms of its scattering properties.


Introduction
The problem of quantum transport is an ubiquitous topic of non-equilibrium statistical physics. Transport phenomena are often discussed in the linear response context, where they essentially reduce via Green-Kubo theory to studies of dynamical correlations in equilibrium. On the other hand, genuine non-equilibrium treatment [1,2] requires a subtle control over the thermodynamic limit or, more precisely, over an infinite number of degrees of freedom composing the so-called thermal (or particle/magnetic) baths (or reservoirs).
In the non-equilibrium protocol, which has become known as the inhomogeneous quench problem, one prepares an initial state of an infinite low dimensional system (say, a one-dimensional chain) which is out of equilibrium and is characterised by some spatial profile of local physical observables, such as density, temperature, magnetisation, etc (see [3,4] for recent reviews and references). Usually, for simplicity, one considers an initial step profile protocol in which the initial state is given as a product of two distinct equilibrium states for the left and the right semi-infinite part of the system with respect to some point of origin, however more general (say, smooth) profiles can also be considered. It has been shown that a quantum quench starting from such types of inhomogeneous initial states allows for an exact derivation of transport properties in closed-form expressions, either in conformal field theory framework [5][6][7][8][9], where it results in a Stefan-Boltzmann radiative transfer, or for free fermions [10][11][12][13][14][15][16][17][18][19], or even in general integrable systems, where the problem can be effectively treated in terms of the so-called Generalised Hydrodynamics [20,21]. The quantum dynamics after such an inhomogeneous quench has been recently observed in ultra-cold atom experiments [22].
The quantities of central interest for the study of quantum transport in such inhomogeneous out-of-equilibrium settings are the particle current and density (or current and density of any other quantity that satisfies a local conservation law) at fixed position and large time, as well as their asymptotic profile at large times and distances with fixed ratio of distance over time. For ballistic transport, which is the typical case for non-interacting or integrable systems, these profiles are characterised by the formation of two (or more, for multiple types of particle excitations) fronts, moving ballistically to the left and right direction, and the relaxation to a Non-Equilibrium Steady State (NESS) in the middle region between the oppositely moving fronts [23]. A NESS is a statistical ensemble that gives the large time values of local observables at fixed position. More generally, it has been shown that a different ensemble can be associated with the asymptotics at each ray of fixed distance/time ratio [15].
As the inhomogeneous quench protocol intrinsically breaks the translational invariance of the problem, it is natural to ask what are the effects of translational symmetry breaking terms in the Hamiltonian, say point (or localised) defects at the contact between the two semi-infinite systems, prepared in two distinct equilibrium states as in the case of the step profile quench protocol. Such defects can be interacting or non-interacting, integrability-preserving (such as in the famous Kondo [24] or Anderson impurity problems [25]) or integrability-breaking [26,27]. The latter could be considered as minimal models of a quantum Kolmogorov-Arnold-Moser theorem scenario. General ideas for the treatment of transport in impurity models have been put forward in the past. At low temperatures Kane and Fisher [28] studied a Luttinger liquid with a weak link and found a non-trivial phase diagram, where RG flows either to a vanishing or a unit quantised conductance. Several approaches have been developed based on integrability [29][30][31] (see also [32] for a review and [33,34] for more recent discussions) and more recently using conformal field theory [35]. A bosonisation treatment has also been used to study energy transport after joining two XXZ spin chains with different parameter values [36]. Another approximate approach is based on Landauer's theory [37] (see also [38] and references therein for related numerical works), which refers to quantum transport through an impurity in contact with two thermal baths at different chemical potential. The Landauer-Büttiker formula describes the steady state current through the impurity as a function of the voltage (Fermi level difference). However, explicit results for the unitary dynamics of the whole closed system far from equilibrium (i.e. for strong inhomogeneity bias) are rare and limited to hard-wall defects [15,39,40].
On the other hand, many analytical studies of the inhomogeneous quench problem are based on the so-called semiclassical or hydrodynamic approach. This approach and its various extensions have been widely used mostly in essentially non-interacting systems, in both homogeneous [41][42][43][44][45][46][47][48] and inhomogeneous settings [13,[17][18][19]49], but more recently also in genuinely interacting integrable systems [20,21,48,[50][51][52][53][54]. In defect related out-of-equilibrium problems it has been applied to the special cases of a quench in the presence of a hard-wall boundary [15] and of a local quench induced by a moving defect [55]. Even though it is not a priori obvious why this approach should be anything more than an approximation, it is conjectured that it gives correct predictions for the exact asymptotics of physical observables at large times and distances, as it turns out to be in perfect agreement with all available numerics.
In this work we study inhomogeneous quenches in the presence of non-interacting defects, which can be studied in an exact and mathematically rigorous way. We focus on a onedimensional lattice system of free hopping fermions (equivalent to the XX model of spin 1 2 ) with a non-interacting particle-number preserving defect. Our goal is to derive from first principles the exact asymptotics of the particle current and density (not only close to the defect, but in the whole system), which we do using a new analytical approach that handles rigorously the thermodynamic and large time/distance limit. In this way we test the validity of the semiclassical approach, demonstrating that it reproduces correctly the exact results away from the defect, but also identifying in certain cases corrections close to the defect that cannot be captured by it. Such corrections are present, for example, when the defect produces localised bound state excitations. These corrections consist in persistent oscillations of physical observables, whose strength decays exponentially with the distance from the defect, and in different time-averaged values in comparison with those away from the defect. As physically expected, the defect obstructs the motion of particles diminishing transport. As a result, we find that the initial density difference on the left and on the right side is partially preserved in the steady state. All of the above can be clearly seen in Fig. 1 showing spacetime plots of the particle density and current for a typical case of such an inhomogeneous quench discussed in detail below. Our analysis shows that the dependence of the current/density profiles on the defect strength can be expressed solely in terms of its scattering properties, i.e. the reflection or transmission probability that characterises it.
The paper is organised as follows: In Sec. 2 we present the semiclassical or hydrodynamic approach to quantum quenches and apply it to the present problem finding simple formulas for the asymptotics of particle density and current. We also discuss the analogy with Landauer's theory, demonstrating the compatibility of the semiclassical results with the Landauer-Büttiker formula in the special case when the initial momentum densities on the left and right side are thermal. In Sec. 3 we proceed to the exact analytical calculation of the asymptotics of current and density profiles. In Sec. 4 we compare our results with numerics finding perfect (numerical data). The initial state has different density on the left and right side and the dynamics is that of a free fermion lattice system in presence of a hopping defect at the origin (numerical data for ratio of defect over bulk hopping j = 1.2; values rescaled between 0 and 1; see Fig. 6 for line plots). Notice the ballistic motion of the fronts, the partial preservation of density difference between the two sides of the defect and the finite distance effects close to it (density corrections and persistent oscillations that decay exponentially with distance from the origin). agreement. In Sec. 5 we generalise the analytical calculation to any non-interacting defect and compare with the semiclassical formulas. We summarise our conclusions in Sec. 6. Effects at finite distance from the defect, partially due to the presence of bound states, are discussed in App. A. Some mathematical tools that are crucial for the rigorous derivation of the results are presented in detail in App. B.1 and B.2.

Semiclassical or Hydrodynamic approach
The asymptotics of local observables after a quantum quench can be successfully captured by a simple quasiparticle picture, known as semiclassical or hydrodynamic approach. In this section we present a description of this approach and apply it to the present problem of an inhomogeneous quench in the presence of a non-interacting defect, deriving explicit expressions for the current and density asymptotics. We later compare these with the exact results and confirm their validity.
According to the semiclassical approach, at least as far as the asymptotic behaviour of the system is concerned, its dynamics can be derived by considering a statistical system of classical quasiparticles that are produced at the time of the quench and subsequently travel ballistically throughout the system. The quasiparticles correspond to those that would describe the excitations of the bulk part of the post-quench Hamiltonian and their velocities are equal to the corresponding group velocities of excitations. In the present problem the bulk Hamiltonian is non-interacting i.e. there are no quasiparticle collisions which could modify or "dress" their velocities.
For inhomogeneous quenches, the initial probability distribution ρ 0 of the created quasiparticles depends on both their position x and momentum k and is deduced for each position from the local values of the quench parameters. More specifically, for a step-like initial state we have where the momentum distributions ρ 0L/R (k) denote the densities that correspond to homogeneous quenches with initial states given by the left and right side restriction of the actual inhomogeneous initial state. To be specific we will focus on a special initial state with flat densities but the following results are valid in general. The time evolution of local observables, like the density n(x, t), is then calculated tracing the classical trajectories of the quasiparticles and integrating the number of those reaching the point x at time t in all possible ways. In absence of the defect, this means integrating over all initial points x from which these quasiparticles may originate. In the presence of the defect on the other hand, we should also take into account the fact that quasiparticles coming from the opposite side have first been transmitted through the defect and some of those coming from the same side can now get reflected off the defect and change direction. These scattering effects amount to introducing weights for the contribution of those groups of particles, equal to the transmission and reflection probabilities T (k) and R(k) respectively. The two cases are schematically represented in Fig. 2.
We first focus on the simpler case where no defect is present and next study the modifications introduced by the presence of the defect.

Absence of defect
According to the above arguments, in the absence of the defect the density is given by where v(k) is the group velocity, which in the non-interacting case is given by the derivative of the dispersion relation v(k) = dE(k)/dk. Performing the integral over the initial quasiparticle Already at this step, we observe that the semiclassical picture captures several of the general characteristics of the asymptotics. First, due to the step-like form of the initial state, the density n(x, t) (and similarly all other local observables) is actually a function of a single variable, the "ray" ratio ξ = x/t, instead of x and t separately. Second, the value of joint positionmomentum density on each of the rays of fixed ratio ξ is given by a ξ-dependent momentum distribution, i.e. with In the simplest case of a monotonously increasing function v(k), ρ ∞ (k; ξ) has a simple stepwise form: for v(k) > ξ it is equal to the initial momentum density of the left side of the system, while for v(k) < ξ it is equal to that of the right side (Fig. 3). The threshold momentum k * separating the two contributions for each ray is therefore the ξ-dependent solution of the equation v(k * ) = ξ. At the middle, ξ = 0, this reduces to the usual observation that the NESS is built out of rightmoving modes carrying the initial density ρ 0L (k) of the left side (where they come from) and left-moving modes carrying the initial density ρ 0R (k) of the right side. The current J(x, t) can be calculated in a similar way. Taking into account that in the semiclassical approach it is given by we finally find Alternatively, we can derive the last result directly from the formula for the density using the continuity equation. Indeed in the scaling limit the continuity equation becomes ∂ t n(x, t) = −∂ x J(x, t) from which we find J(x, t) = −´x −∞ ds ∂ t n(s, t) and due to the scaling form of the density Substituting (3) and (4) we obtain (6).
In the present problem of the XX spin 1 2 model, the group velocity is This is not a monotonous function in the interval [−π, +π], however it is monotonous in [−π/2, +π/2] and has symmetric values outside of the interval. Therefore, as long as the initial densities ρ 0L/R (k) are even and symmetric with respect to ±π/2, we can restrict the momentum integral in (3) between −π/2 and +π/2 and write The solution of (5) in this restricted interval for v(k) given by (8), is Using this solution, we find In the special case of constant initial densities (2) we obtain

Presence of defect
We now focus on the case where a defect is present. In the semiclassical approach the quasiparticles are still the same free particles and have the same initial densities as before, but now they can also scatter with the defect. This means that when calculating the time evolution of the density n(x, t), there are two modifications with respect to the previous case. First, the number of quasiparticles with momentum k that come from the opposite side is reduced by a factor equal to the transmission probability T (k), since they have been transmitted through the defect. And second, there are now quasiparticles originating from the same side that have been reflected off the defect with probability R(k) and changed direction of motion. We therefore find that the time evolution of the density on the left side is given by In the above we used the fact that transmitted quasiparticles keep moving in the same direction with the same velocity, while reflected ones change direction. Substituting (1) and using the fact that ρ 0L/R (k) are even function we obtain Working similarly for the right side, we find Therefore there is a difference between the asymptotic values on the left and on the right of the defect, which is We similarly find that the asymptotic profile of the current is given by which is valid on both sides, as we can easily verify using the relation R(k) + T (k) = 1 and the fact that ρ 0L/R (k) and T (k) are even functions while v(k) is odd. In the same way we can check that the above expression is symmetric under x → −x.
The current at the middle is In the special case of constant initial densities (2) we find that the current at the middle is Note that the current is proportional to the difference of the initial density on the left and right halves and has a direction from the side of higher to that of lower density. Also note that the current vanishes when the defect has infinite strength, blocking completely the passage of particles from one side to the other (T (k) = 0) and takes its maximum value when no defect is present (T (k) = 1). This is of course what we expect from physical intuition, since the defect tends to obstruct the motion of particles and therefore to diminish the value of the current and to partially preserve the density difference between the left and right side that is inherited from the inhomogeneous initial state.
For the hopping defect which we are going to study in more detail, the reflection and transmission probabilities turn out to be where j is the ratio of the defect hopping over the bulk hopping and measures the defect strength. These expressions can be readily found by solving the quantum mechanical problem of an incident particle travelling with momentum k towards the defect and being reflected and transmitted through it. Notice that these expressions are clearly symmetric under inversion of the defect strength j → 1/ j. According to the above semiclassical analysis, this symmetry will be manifest in all results describing the asymptotics of density and current.

Comparison with Landauer's theory
We will now compare the inhomogeneous quench problem with that of an open and typically finite quantum system coupled to two thermal baths at different chemical potential, which was studied by Landauer [56]. This theory describes the current through the system as a function of the voltage, i.e. chemical potential imbalance, and is relevant for the physics of quantum dots and quantum wire junctions. This problem is different from the quench in terms of physical settings and questions: In the quench problem the system is extended and closed (that is, the dynamics is unitary) and we are interested in the spatial profile of observables which are timedependent for all times. In Landauer's problem the system is open, the internal dynamics of the baths is ignored, as they are assumed to have fixed thermal spectral densities, and we are interested in the time-independent conductance as a function of the voltage, the Fermi level difference of the external baths. Despite these differences, there is a strong physical analogy between the two problems [37], since one may interpret the quench problem as follows: The defect can be viewed as an open system coupled to the two halves of the rest of the system. At large times, assuming that a NESS is reached, the particle current coming towards the defect from each side is constant and determined by the particle density produced at quench time far from the defect on the left and right side. Therefore the two halves can be viewed as infinite baths with fixed densities in contact with the defect. It is then reasonable to expect that Landauer's theory is compatible with the semiclassical predictions for the NESS current at the middle in the quench problem. An important difference however is that the effective baths in the quench protocol are not thermal, but described instead by Generalised Gibbs Ensembles corresponding to the homogeneous quenches on the left and right side asymptotically far from the origin. This means in particular that the left/right density imbalance is generally non-zero at all energy levels, in contrast to the case of low-temperature thermal states that is usually considered in Landauer's problem. For the same reason the density imbalance cannot be generally described in terms of voltage, i.e. difference between Fermi levels of thermal distributions.
It is easy to see that our semiclassical expression (13) for the current at the middle is compatible with Landauer's results in the special case of thermal initial densities with different Fermi levels on the two sides. Indeed, changing the integration variable from momentum to energy using v(k) = dE/dk, taking into account the double energy level degeneracy and finally which is the well-known Landauer-Büttiker formula. At arbitrary voltage but zero temperature 1/β → 0 we obtain which is another well-known result of Landauer's theory.

The quench protocol
We consider a system of free fermions in a one-dimensional lattice of length L, assumed for simplicity to be even. We prepare the system in an initial state ρ 0 that is diagonal in the local basis with a stepwise inhomogeneity in the density and we let it evolve under the free hopping fermion Hamiltonian with a defect in the hopping located at the origin (precisely, between the sites 0 and 1) As can be seen, we have set the bulk hopping J b = 1 and the defect hopping J d = j.
We are interested in the expectation value of the current, and the particle density, in two types of limits: i) at fixed position x and large time t → ∞, and ii) at large times and distances, but keeping the ratio x/t fixed. More generally we are interested in the asymptotics of the two-point correlation function from which the above and any other observable can be derived. This is because both the initial state and the time evolution are Gaussian in the fermionic modes c x , c † x , therefore the problem reduces to a single-particle problem and multi-particle observables can be derived using Wick's theorem.

Post-quench eigenstates
We start by characterising the eigenfunctions φ E (x) and corresponding eigenvalues E of the Hamiltonian H in the single particle sector. These are either odd or even with respect to the spatial reflection x → 1 − x, so we can distinguish them with a sign index σ = ± To determine their form, we can focus on their restriction in the right semiaxis (x ≥ 1) and use the condition that they are eigenvectors of the single particle Hamiltonian together with reflection symmetry in order to find the allowed discrete values of the energy E.
In the bulk (anywhere except at the origin) the eigenvector equation becomes This is a recursive equation with general solution where the new "momentum" parameter k is given as a function of E through which is the dispersion relation. The parameter δ σ (k) is a phase-shift or "extrapolation length" to be determined by the "initial" conditions of the recursive equation, that is, by the conditions at the origin. From now on we will be using k instead of E to parametrise the eigenfunctions as φ σ (k; x). At the origin the eigenvector equation becomes which, together with the reflection symmetry (22) that can be written as a condition on the middle site values φ σ (k; 0) = σφ σ (k; 1), gives This plays the role of "initial" condition for the above recursive equation (23). By requiring that the general solution (24) satisfies it, we find that δ σ (k) must satisfy the equation which defines δ σ (k). At the right boundary x = L/2, the eigenfunctions satisfy the equation which can be written as cos k (L/2 + 1 − δ σ (k)) = 0, and constitutes the quantisation condition determining the discrete values of the parameter k (or E, through (25)). Lastly, the normalisation factor N turns out to be For j > 1 there are also two bound states, one even and one odd, that are localised at the defect and correspond to imaginary k. In the thermodynamic limit L → ∞, they are given by with κ = log j and the corresponding energies

Calculation of the two point correlation function
We will now calculate the two point correlation function in the thermodynamic limit and derive its asymptotics at large times, using a mathematically rigorous approach [57]. Let us focus on the case j < 1, in which there are no bound states. The calculation in the case j > 1 is analogous and is outlined in App. A. The time evolution of the two-point correlation function can be expressed in terms of the eigenfunctions as where is the free fermion propagator, and C 0 ( y, y ) = Tr ρ 0 c † y c y is the initial correlation function. For our initial state (19) with only diagonal correlations in the local basis, this is given by C 0 ( y, y ) = δ y, y n 0 ( y), and using the stepwise form of the initial density (19), we find The first term is irrelevant for transport, since by straightforward calculation using the completeness of the eigenstates we find Using the reflection symmetry of the eigenstates, the second term can be written as We calculate the spatial sum of the last expression and eliminate L using the quantisation condition (27) in the form e ik(L/2+1−δ σ (k)) = −e −ik(L/2+1−δ σ (k)) in order to perform the thermodynamic limit L → ∞ easier Having eliminated L, the summand of (35) is now a function of variables k, k that, in the thermodynamic limit, spread uniformly in the interval from 0 to π, therefore we can write the sums over k, k as integrals. However, since the functions δ σ (k), δ −σ (k ) are implicitly known and defined in terms of the discrete values of k, k , we do not know if they converge uniformly to smooth functions in this limit. Luckily, this problem can be circumvented easily: if we multiply with φ σ (k; x)φ −σ (k ; x ) and use the defining equation (26) of these implicitly known functions, we can eliminate them completely Notice that, while the original summand of (35) is well-defined for all discrete values of k, k and any finite L, the continuous function to which it converges in the thermodynamic limit exhibits a pole at k = k . This pole, which reflects precisely the inhomogeneity of the initial state, is going to play a crucial role in the derivation of the NESS, as it turns out to completely determine the combined thermodynamic and large distance and time asymptotics. Its presence prevents us from writing the thermodynamic limit of the sums directly as integrals, since we need to find the correct contour prescription for how to pass around the pole. This can be done using a standard complex analysis trick, described in more detail in Appendix B.1: we introduce a meromorphic function with poles of unit residue at the discrete values of k that we sum over, so that we can write the sum over k as a contour integral along straight lines just above and below the real axis, but paying attention to subtract the residue of the additional pole at k = k. Then we can safely take the thermodynamic limit, in which only the contribution of the extra pole and of one of the above straight line integrals survives.
Using this trick (equation (64) of the appendix) and substituting into (34), we finally obtain an exact formula for the thermodynamic limit of the correlation function where the integral over k runs along a contour from 0 to π slightly below the real axis. The last expression holds only for x, x > 0, since in the above calculation we used the right half side expression of the eigenfunctions φ σ (k; x). The explicit expression in different spatial regions can be deduced from the above, as already mentioned, by exploiting the reflection symmetry with respect to the origin.

Large time fixed position density and current
In the large time limit t → ∞ at fixed positions x, x the contour integral over k in (37) vanishes, since from the dispersion relation (25) we find that the exponent of e −iE(k )t has a negative real part −(dE/dk )εt = − sin k εt < 0 for all k from 0 to π. Therefore only the residue of the pole at k = k survives and we find (valid for x, x > 0). (38) Notice that the multiple sums over eigenstates and the sum over all lattice sites in the original expression (34) have been reduced to a single integral over momenta in the combined thermodynamic and large time limit.
The last expression (38) for the correlation function determines completely the NESS that forms at large times. In particular we find the large time value of the current at the origin as a function of the defect strength For j > 1, as shown in App. A, the large time current is given by the above expression plus additive corrections which are due to the presence of bound state excitations. We see that the bound states induce persistent oscillations in time at frequency 2( j + 1/ j) equal to their energy and which decay exponentially with the distance from the origin at a characteristic length 1/(2 log j).
Notice the symmetry of the expression (39) under the discrete transformation j → 1/ j. This reveals an interesting strong-weak duality which, as we will see in the following, is manifest in all results describing the asymptotic behaviour of the system and which is neither present nor obvious in the microscopic description of the problem. As anticipated in Sec. 2, the reason for the emergence of this symmetry is that the asymptotic behaviour can be fully determined from the quasiparticle reflection or transmission probabilities that characterise the defect and which for the hopping defect turn out to be symmetric under inversion of the defect strength.
Another observation we can derive from the expression (38) is that for j < 1 the current J(x, t) at fixed position and large time is independent of the position lim t→∞ J(x, t) = lim t→∞ J(0, t) (valid for j < 1). This is because the expression for lim t→∞ J(x, t) is given by the same integral as lim t→∞ J(0, t) independently of the position x.
The large time density n(x, t), on the other hand, exhibits a stepwise inhomogeneity in space but is otherwise independent of the position. More specifically, we find that for x ≥ 1 while for x ≤ 0, using the reflection symmetry we have The calculation of the above integral is a standard application of residue theorem: we first express it in terms of the complex variable z = e 2ik , separate the integrand in powers of z x and note that the various terms have poles at either z = j 2 or z = 1/ j 2 or both, which are located inside or outside of the unit circle depending on whether j < 1 or j > 1. In the present case j < 1, it turns out that these poles do not contribute to the above expression. After some algebra, we finally find that for j < 1 the NESS density exhibits a simple step but is otherwise independent of the position For j > 1, as shown again in App. A, the density is given by the above expression plus additive corrections i.e. apart from a steplike inhomogeneity at large distances, there are also corrections that decay exponentially with the distance from the origin together with persistent oscillations in time that also decay away from the origin. In this case the large time limit does not exist strictly speaking, due to the presence of the oscillatory terms. The density difference on the left and on the right of the defect and far from it, is given by We observe that the asymptotic NESS density difference is proportional to and has the same sign as the initial density difference ∆n 0 = n 0 (x > 0) − n 0 (x < 0) = 2µ, that is the defect preserves part of the initial density difference. Note that, even though for j < 1 the difference on the two sides far from the origin is equal to the difference between the sites 0 and 1, this is not true for j > 1 where, as already mentioned, there are corrections close to the origin that result in a time-averaged density difference equal to where the bars denote infinite time averages. Note that this difference has the opposite sign of the initial density step.

Large distance and time asymptotics of density and current
Apart from the asymptotics at fixed position x and large times t, we can also derive the asymptotics at large times and distances from the origin keeping their ratio ξ = x/t fixed. To this end we go back to the expression (37) for the correlation function in the thermodynamic limit, expand the integrand in terms containing exponentials of the form e ±ik x−iE(k )t and derive the asymptotics of the k -integral of each term separately. This can be done by deforming if necessary the integration contour slightly above the real k -axis to make sure that it passes only through regions where the exponent has a negative real part and therefore the new integral decays in the considered limit. The necessary contour deformation depends on the value of x/t. If the pole of the integrand at k = k is crossed during the contour deformation, then we have to take into account the contribution of its residue. This complex analysis method is described in more detail in Appendix B.2. Using this method (equation (65) of the appendix), we finally find that the density n as a function of the rescaled variable x/t is where F ( j; ξ) is the following piecewise defined function Note that the above function satisfies the expected properties: It equals ±1 outside of the lightcone |ξ| = |x|/t > v max = 2 on the left/right side respectively, reducing to the initial left/right half side density values. For j = 1 i.e. in the absence of defect, it reduces to which describes the rescaled density profile for an inhomogeneous initial state evolving under the homogeneous defectless Hamiltonian, derived in [10]. For j = 0 or ∞, i.e. when the defect is impenetrable and keeps the two halves of the system completely disconnected from each other, F becomes that is, there is no transport at all and the density remains everywhere equal to the initial value.
For any j = 1 on the contrary, it exhibits a discontinuity at ξ = 0 which is the most characteristic effect of the presence of the defect. Note that the value of the discontinuity jump is equal to the value ∆n NESS in (44), showing that these two different limits match. The asymptotic profile of the current J as a function of the rescaled variable x/t can be calculated directly from the density profile n(x/t) using the continuity equation where the function G( j; ξ) is defined as For j = 1 we recover the known result for the case of absence of the defect 4 − ξ 2 for − 2 < ξ < 2 0 otherwise, and for j = 0 or ∞, we have G(0, ξ) = G(∞, ξ) = 0, as expected. Note also that the current is continuous as |x|/t → 0 and its value is given by that gives the same value (39) found earlier for the large time current at the origin J NESS , which means that the two limits give the same result, as for the density profile. Lastly note that the functions F ( j; ξ) and G( j; ξ) are symmetric under inversion of the defect strength j → 1/ j, i.e. F ( j; ξ) = F (1/ j; ξ) and G( j; ξ) = G (1/ j; ξ). This shows that the strong/weak defect duality pointed out earlier for the NESS expectation values (4) and (44) is manifest also in the large distance and time asymptotics of the density and consequently the current.

Comparison with numerics
We will now compare our analytical results with numerical data. Since the defect considered here is non-interacting, the full quench dynamics can be obtained efficiently by exact diagonalisation 1 . The large time and distance asymptotics can be observed by considering a sufficiently large finite system and long times up to some order smaller than the system size divided by the maximum group velocity, in order to exclude finite size effects.
In Fig. 4 the analytical results (39) and (44) for the NESS current and density difference are plotted as functions of the defect strength j together with numerics. The numerical data are extracted from the current and density profiles at the largest time reached in the computation (t = 50). For j > 1, in order to subtract the time oscillation corrections (40) and (43) close to the defect, we average over a sufficiently large time window (from t = 40 to 50) to integrate them out. For the estimation of the NESS density difference we also apply a linear fit on the spatial profiles in order to extrapolate the values on the left and right from the origin, eliminating as much as possible the exponentially decaying corrections.
There are some small deviations between analytics and numerics for the density difference only for j > 1 and close to 1. These are exactly due to the above finite distance effects and, more precisely, due the fact that they spread over a characteristic length that diverges as j → 1 + , meaning that at any finite time there is some mixing between these features and the asymptotic form of the profile. For this reason, data extrapolation from linear fitting is less effective in this region and the estimates of ∆n NESS are distorted towards ∆n d given in (45) which has opposite sign. Except for these clearly understood deviations, we observe perfect agreement between analytical and numerical results. Note also the reflection symmetry of the curve with respect to the vertical line at j = 1 in logarithmic-linear plot, which demonstrates the emergent strong-weak duality symmetry j → 1/ j.
The analytical expressions (46) and (48) are plotted against numerical data in Fig. 5 and  6 for several values of j ≤ 1 and j ≥ 1 respectively. The numerical data correspond to the maximum time reached in our computation (again t = 50).
Notice the oscillatory behaviour of the numerics close to the fronts. This is shown in more detail in Fig. 7 where numerical data for the current profile close to the right front are plotted . Current and density units are −4µ/π and 2µ respectively, with 2µ the size of the initial density step. Analytical results (red curves) given by (39) and (44) vs. estimates extracted from numerics (blue dots).
for several different times together with the analytical curve. These features of the profile are described by the Airy function [16,18] and are finite time corrections, since even though they broaden in space as time increases, their width scales slower than linearly with time and therefore shrinks to a point as t → ∞ when plotting in terms of the rescaled variable x/t. This is the reason why this effect is not present in the asymptotic expressions (46) and (48). Also notice the oscillatory and exponentially damped deviations of the numerics close to the defect for j > 1 (Fig. 6), which, as already explained, are restricted to finite distances from the origin and therefore also disappear in the rescaled variable plots when t → ∞.
With the exception of these expected corrections due to finite distance or time effects, the analytical curves match perfectly with the numerics. In Fig. 8 later on we perform the same comparison between analytics and numerics for a different type of non-interacting defect, observing once again perfect agreement.

General defect: asymptotics from the scattering matrix
It is easy to verify that the above analytical results (46) and (48) for the asymptotic density and current profiles are identical to the semiclassical predictions (9), (10) and (12) once the reflection and transmission probabilities are substituted by their explicit form (15) and (16) for a hopping defect. We will now show that this is generally true for any type of non-interacting and particle number preserving defects. To this end we are going to generalise our analytical calculation using the scattering matrix that completely characterises such a defect. Scattering matrix based analytical approaches for the study of quantum transport in presence of a defect have been used in the context of quantum wires or quantum dots coupled to infinite baths [58][59][60][61] and recently for the study of light cone effects after a local quantum quench [62].  (46) and (48).
We consider a general non-interacting defect V d , located in a finite spatial interval [− /2 + 1, + /2] centred at the origin (here x = 1/2), focusing for simplicity on reflection symmetric defects The eigenstates are plane waves with a phase shift between the left and right side of the defect, which can be written as The scattering or S-matrix is defined through the relation and r L/R , t L/R are the reflection and transmission probability amplitudes for a particle incident from the left (D = 0) or right (A = 0) respectively. For defects symmetric under spatial reflections, the S-matrix is a symmetric matrix with only two complex parameters From the requirement of particle conservation, the S-matrix must always be a unitary matrix, i.e. S † S = I, which means that r, t satisfy the relations r t * + r * t = 0.
As always the reflection and transmission probabilities are defined as the norm squares of the corresponding amplitudes R = |r| 2 , Note also the useful formula that can be found from the above. The eigenvalues of the S-matrix are the unitary numbers t ± r.
In the symmetric defect case, the eigenfunctions are either odd or even We therefore find additional relations between the coefficients A, B, C, D We can choose to parametrise the solutions (49) with respect to the right side coefficients C, D. Moreover it is convenient to rewrite them in trigonometric form as = e ik (r * (k) + σt * (k)) .
The quantisation conditions can be still written in terms of δ σ (k) as Having expressed the eigenstates in the same form as in the previous special case, the calculation of the correlation function and the asymptotic profiles can be done by repeating the method of Sec. 3.3 and 3.4. We thus find which can be equivalently written in terms of the phase shifts e 2ikδ ± (k) as From these expressions we can derive formulas for the NESS current and density difference far on the left and right of the defect, taking into account that potential presence of bound states or singularities of the S-matrix would only result in corrections that decay exponentially with the distance from the defect and therefore vanish far from it. Using (52), it turns out that both of these values can be expressed in a very simple form in terms of only the reflection R(k) or transmission probability T (k) = 1 − R(k) Lastly we can derive the density and current asymptotics at large times and distances, which are also expressible in terms of T (k) only As already mentioned, we can verify that all above expressions reduce to the previously derived formulas for the special case of hopping defect. The reflection and transmission probability amplitudes are resulting in the corresponding probabilities (15) and (16). Substituting r(k), t(k) from (59) and (60) we recover all previously found formulas for the asymptotics of the correlation function, density and current.
As an additional verification test we now consider a different type of non-interacting defect, a density defect at the sites 0 and 1 with λ the strength of the defect. We can easily derive the reflection and transmission coefficients  (57) and (58) with R(k) given by (61).
resulting in a reflection probability Fig. 8 shows the comparison between analytical and numerical profiles for various values of the defect strength λ. Only positive values are used since the asymptotic profiles for values of λ with opposite signs turn out to be identical. This is due to the symmetry of R(k) in (61) under the transformation (λ, k) → (−λ, π − k). The agreement between analytics and numerics is once again perfect everywhere except close to the origin, where as expected finite distance effects due to the defect cause deviations of finite time numerical results from the asymptotic values. As a last note, we should comment on the applicability of our analytical method for general inhomogeneous initial states. Even though we have focused our analysis on initial states of the product-state form with a sharp density step between the two adjacent middle sites as described by (19), our method applies to any other initial state that is characterised by different asymptotics on the left and right side far from the middle (independently of its behaviour in the middle) and that is Gaussian with respect to the fermionic fields. This class of states includes not only product-states with smooth density profiles n 0 (x), but also non-product states described by two-point correlations C 0 (x, y) such that lim r→∞ C 0 (x − r, y − r) = C 0L (x − y) and lim r→∞ C 0 (x + r, y + r) = C 0R (x − y) with C 0L (x) = C 0R (x). To see this, one has to look at our derivation in more detail: What we essentially showed is that the large time and distance asymptotics of correlations is completely determined by the location and residue of the momentum-space pole of the correlation function (37) at equal momenta. The origin of this pole can be traced back to the step-like form of the initial state and its residue can be shown to be determined by only the large distance asymptotics of the initial state correlations and not by the slope of the step or any other intermediate spatial features. The only other relevant assumption in our method is that the momentum-space expression for the correlation function is free from singularities on the real momentum axis, which is guaranteed at least for initial states with exponentially decaying correlations.
On the other hand, typical examples of inhomogeneous initial states in spin chains may not necessarily fall in the category of Gaussian fermionic states (due to the non-linearity of the Jordan-Wigner transformation from spin to fermion variables). For non-Gaussian initial states, the two-point correlation function is no longer sufficient to provide a complete description of the state. However even in this more general case the large time asymptotics under a free fermion Hamiltonian is expected to be described by a Gaussian steady state that keeps only memory of the initial two-point correlations (provided that initial correlations satisfy the clustering property [63,64]), therefore the same arguments apply.

Conclusion
In this work we studied the effects of a non-interacting defect on quantum transport after inhomogeneous quantum quenches. We showed in an exact and rigorous way that the asymptotics of particle density and current are correctly given by a simple semiclassical approach, except for finite distance corrections close to the defect. Moreover, we derived exact expressions for these finite distance effects that cannot be captured by the semiclassical approach. Knowledge of their general characteristics (time oscillations and exponential decay with distance) helps to distinguish them from the physically more interesting behaviour of the large distance asymptotics.
An obvious open question that arises is how to handle genuinely interacting defects and whether the previous semiclassical approach can be suitably generalised to describe the asymptotics in this more general case. Our analysis of the non-interacting case suggests that one prerequisite for the exact solution of the interacting problem is to understand the scattering effects induced by the defect on multi-particle states. We leave this exciting problem for future work.
Lastly, it is worth emphasising the analogy between an inhomogeneous quench and the Landauer problem. In both problems the physical interpretation and the mathematical treatment are similar. This analogy points to a broader connection between out-of-equilibrium physics of open and closed extended systems that would be interesting to formulate more rigorously and possibly exploit in order to study more complicated problems in both contexts.  45)) between the two defect sites 0 and 1 that is inverted in comparison with the asymptotic averaged difference ∆n NESS (Eq. (44)) far from the defect. In contrast, for j < 1 the jump between the defect sites is exactly the same as that far from the defect and the density profile is independent of the position, apart from this jump.

B Complex analysis tricks B.1 Conversion of sums into integrals in the presence of a pole singularity
We consider a sum of the form b k=a step: δk where f (k) is regular at all points k n = a + nδk, n = 0, 1, 2, ... but may have singularities at other points in the interval (a, b) and we want to write it as an integral in the limit δk → 0.
We first write the sum as a sum of contour integrals on infinitesimal circles encircling each of the numbers k n with a suitable integrand 1/(e 2πi(k−a)/δk − 1)/δk that has at those points first order poles all with residue 1/2πi. Then we merge the contours into a single contour enclosing the straight line interval in which all numbers k lie (i. e. [a, b]), but making sure to remove any poles of the function f (k) in this line interval by subtracting their residues. Then we take the δk → 0 limit, in which one of the two lines of the contour, the one below, gives a vanishing contribution since 1/(e 2πi(k−a−iε)/δk − 1) → 0. The one above instead gives simply an integral along the real line interval [a, b] shifted by an infinitesimal imaginary number +iε, since 1/(e 2πi(k−a+iε)/δk − 1) → −1.
We therefore find that the original sum converges to an integral along a contour just below Figure 11: Complex analysis method for the derivation of the asymptotics of the integral.
The condition on which side of the real axis the contour should be shifted at some value of k is Re [i((k − iε)ξ − E(k − iε))] < 0, i.e.
Therefore the points at which the shift of the contour should change sign are given by the values k * that are solutions of the equation E (k * ) = ξ.
If ξ > E (q) then the pole at k = q is inside a region where the direction of the deformation should be switched to ensure exponential decay of the integral and therefore the pole does contribute, otherwise not. We finally find that b−iε a−iε dk g(k)e +i(kξ−E(k))t t→∞ −−−−−→ x/t: fixed 2πi e +i(qξ−E(q))t Θ(ξ − E (q))Res k=q g(k).
Assuming, for example, that E (k) is a monotonously decreasing function, there is only one solution k * (ξ) and the pole contributes only when q > k * (ξ) (Fig. 11).