Extensive Long-Range Entanglement in a Nonequilibrium Steady State

Entanglement measures constitute powerful tools in the quantitative description of quantum many-body systems out of equilibrium. We study entanglement in the current-carrying steady state of a paradigmatic one-dimensional model of noninteracting fermions at zero temperature in the presence of a scatterer. We show that disjoint intervals located on opposite sides of the scatterer, and within similar distances from it, maintain volume-law entanglement regardless of their separation, as measured by their fermionic negativity and coherent information. The mutual information of the intervals, which quantifies the total correlations between them, follows a similar scaling. Interestingly, this scaling entails in particular that if the position of one of the intervals is kept fixed, then the correlation measures depend non-monotonically on the distance between the intervals. By deriving exact expressions for the extensive terms of these quantities, we prove their simple functional dependence on the scattering probabilities, and demonstrate that the strong long-range entanglement is generated by the coherence between the transmitted and reflected parts of propagating particles within the bias-voltage window. The generality and simplicity of the model suggest that this behavior should characterize a large class of nonequilibrium steady states.

Within the broad field of quantum many-body physics, the study of nonequilibrium phenomena is becoming increasingly intertwined with the analysis of entanglement witnesses.In particular, the scaling of various entanglement measures with the size of a subsystem quantitatively captures canonical nonequilibrium behaviors, such as thermalization [1][2][3] or the violation thereof [4][5][6], in closed systems subjected to an initial quench.In quench problems of this type, transient effects of long-range entanglement are signatures of integrability [7][8][9][10], and the dynamics as well as the stationary values of the entanglement entropy, negativity, and mutual information are used for the classification of out-of-equilibrium models and their phases [11][12][13][14][15][16][17][18][19][20][21][22].
This success motivates the examination of entanglement properties also in open systems, and specifically those of their steady states, which may give rise to unique entanglement structures [23][24][25][26][27][28][29].Current-carrying states of inhomogeneous systems offer a promising ground for such an analysis, as recent studies have revealed that they can naturally sustain long-range quantum coherent correlations [30][31][32].In this context, scaling laws of steadystate entanglement measures were shown to be closely related to the localized-diffusive phase transition of the open noninteracting Anderson model [33].In this work we show that nonequilibrium conditions may lead to an even more striking behavior of quantum information measures.This is achieved through the study of an elementary model for an inhomogeneous system in a current-carrying state, where the mechanism underlying its unusual entanglement properties is exceptionally transparent.
Beyond the role of entanglement measures as fundamental quantities, their estimation is inextricably linked to the construction of useful tensor-network simulations of condensed matter systems [34,35].Strong (volume-law) entanglement, which is commonly found in nonequilibrium quantum many-body states [2,13,36], impedes the utility of these simulation methods [37].One possible key for their improvement is thus the uncovering of nontrivial entanglement structures in certain classes of states, like the one that is the subject of this work.Steady states that are predicted to give rise to strong entanglement are also of potential interest from a technological standpoint, as entanglement is an essential resource for quantum information applications [38][39][40].
In this work, we examine the long-range entanglement induced by a current-conserving scatterer in the voltagebiased steady state of a 1D noninteracting fermion system.We treat this problem generally, without imposing specific assumptions regarding the structure of the scatterer other than it being smaller compared to all other length scales.We study the correlations between two disjoint subsystems located on opposite sides of the scatterer: A L on its left, and A R on its right.The total amount of correlations is regularly quantified using the mutual information (MI) between the two subsystems, ( Here A = A L ∪ A R , and S X = −Tr [ρ X ln ρ X ] is the von Neumann entanglement entropy of a subsystem X [41], with ρ X being the reduced density matrix of X.
Given that A is in a mixed state, however, the MI has limitations as a measure of entanglement, since it takes into account both classical and quantum correlations [42].Therefore, we also address the fermionic negativity [43][44][45] between A L and A R , an entanglement monotone defined as where ρ A is obtained from ρ A by applying a partial time-reversal to either A L or A R .Interestingly, our analysis shows that the MI and negativity follow a similar scaling, a scaling which to the best of our knowledge has not been previously observed in a natural physical scenario.
As our main result, we find that both the MI and the negativity scale linearly with ℓ mirror , the number of sites in A L that, under reflection with respect to the position of the scatterer, overlap with sites in A R (see Fig. 1(a) for an illustration).Remarkably, this steady-state extensive entanglement is long-ranged, as the observed volume-law scaling does not decay with the (similar) distance of the mirroring sites from the scatterer.Moreover, the entanglement depends non-monotonically on the distance of either A L or A R from the scatterer.We analytically derive exact formulas for the asymptotic scaling of the MI and the negativity (Eqs.( 8)-( 9)).Additionally, we demonstrate that the coherent information (CI) [46,47], is not only positive (which is impossible classically) when ℓ mirror is large enough, but also grows with ℓ mirror according to a volume law (Eq.( 10)) 1 .The CI is a lower bound to the squashed entanglement [48,49], another rigorous entanglement measure with axiomatically desirable properties [41], which therefore obeys an extensive scaling as well in regimes where I(A L ⟩A R ) > 0. A simple intuitive explanation for these results is provided by considering that the scatterer splits each incoming single-particle wavepacket into two coherently correlated counter-propagating parts (Fig. 1(b)).Each such split wavepacket with energy within the voltage window generates entanglement, since detecting the particle in one subsystem prohibits its presence in the other.As the number of such wavepackets is proportional to ℓ mirror and independent of the distance between the subsystems (Fig. 1(c)), the correlation measures exhibit a similar behavior.
The paper is organized as follows.In Sec.II we introduce the model for the system and its nonequilibrium steady state that are the subject of this work.In Sec.III we report our analytical results for correlation measures in the steady state.We point out the salient features of these results, and support them through comparisons to numerical results (computed for a specific choice of the scatterer).Sec.IV outlines the derivation of the analytical results, and is limited to the conceptually crucial steps in the derivation, while the technical aspects of the process are mostly discussed in the appendices.In Sec.V we conclude and mention potential future directions arising from this work.
Additionally, the paper includes four technical appendices.In Appendix A we derive the two-point correlation function, which served as the basic ingredient in all of our calculations.Appendix B presents the technical details of our computation method for subsystem entropies, from which (as explained in Sec.IV) the asymptotics of the MI and CI can be immediately derived.Appendix C summarizes the derivation of the formula for the fermionic negativity, which is based on the same method.Finally, in Appendix D we complement the numerical results included in Sec.III with additional numerical tests corroborating our analysis.

II. NONEQUILIBRIUM MODEL
We consider a 1D lattice, occupied by noninteracting fermions and connected at its ends to two reservoirs with different chemical potentials, µ L ̸ = µ R , at zero temperature.The lattice is assumed to be of infinite length, and it is modeled as a tight-binding chain that is homogeneous everywhere, except for a small region at the center of the chain, which we dub the scattering region.The Hamiltonian is thus of the form Here c m is a fermionic annihilation operator for the mth lattice site, η > 0 is a hopping amplitude, m = ±m 0 designate the boundaries of the scattering region, and H scat pertains only to sites with |m| ≤ m 0 and breaks the homogeneity, e.g., through modified hopping terms, on-site energies, or side-attached sites.
The scattering region can be associated with a 2×2 unitary scattering matrix [50], defined for any lattice momentum 0 < k < π: The diagonal (off-diagonal) entries of this matrix stand for reflection (transmission) amplitudes; the left (right) column contains the scattering amplitudes for a particle originating in the left (right) reservoir with momentum k > 0 (−k < 0).The squared moduli of the entries correspond to the transmission and reflection probabilities, respectively T (|k|) and R(|k|) = 1 − T (|k|), for a particle originating in either reservoir with momentum k.These scattering probabilities are the sole property of the scatterer on which our analytical results depend.The single-particle eigenbasis of the Hamiltonian is comprised of extended scattering states with energies ε = −2η cos k, and of bound states localized near the scattering region [50,51]; we ignore the latter in our analysis, as they contribute negligibly to correlations between sites that are distant from the scatterer.The current-carrying many-body steady state is pure, with single-particle scattering states originating in the left (right) reservoir occupied up to a Fermi momentum k F,L > 0 (−k F,R < 0), as shown schematically in Fig. 1(a).The Fermi momenta are related to the chemical potentials through Correlation and entanglement measures are calculated with respect to two subsystems A L and A R , each comprised of contiguous sites, with lengths ℓ i and distances d i ≥ 0 (i = L, R) from the scattering region (all of which are assumed to be much larger than the size of the scattering region, 2m 0 + 1): and we also define ∆ℓ i = ℓ i − ℓ mirror .

III. ASYMPTOTICS OF CORRELATION MEASURES
The leading behaviors of the MI and the negativity can be encapsulated by that of the Rényi MI, defined as where n ] are Rényi entropies (which converge to S X as n → 1).We report that, for the nonequilibrium steady state described above, the Rényi MI follows a volume-law scaling with ℓ mirror , where k − = min{k F,L , k F,R } and k + = max{k F,L , k F,R } are the two Fermi momenta that bound the voltage window.The ellipsis (which will be henceforth omitted) represents subleading terms, the dominant of which are logarithmic in the different length scales (ℓ i , d i , and combinations thereof).Full exact expressions for these logarithmic terms can be obtained in the long-range limit d i /ℓ i → ∞ (with d L − d R kept fixed) using methods related to the asymptotic calculation of Toeplitz determinants [52][53][54][55]; the results for these subleading corrections will be discussed in a separate publication [56].
The MI is related to the Rényi MI simply by its definition, through the equality A L :A R , yielding the following asymptotics: The negativity, on the other hand, is not a priori directly related to the Rényi MI, yet our analysis shows that, at the leading (linear) order, The equality I (1/2) is known to arise in quenches of integrable systems [10].Eqs. ( 8) and ( 9) state that, for a generic non-trivial scatterer (i.e., unless ), the MI and negativity both exhibit extensive scaling with ℓ mirror .Additionally, we have found that the CI scales at the leading order as and so it grows linearly with ℓ mirror if ∆ℓ L is fixed.Crucially, the asymptotics in Eqs. ( 7)- (10) do not depend on the magnitudes of d i , and they hold even when d i ≫ ℓ i .That is, the extensive entanglement is long-ranged, and it holds even for subsystems that are very distant relative to their lengths, but that still share mirroring sites.Eqs. ( 8)-( 10) are the central results of this work.The special symmetric case where ℓ L = ℓ R = ℓ and the subsystems are positioned symmetrically relative to the scatterer (d L = d R ) is particularly illuminating with regard to the nature of the steady-state entanglement.In this case we have found that S A scales sublinearly with ℓ, i.e. lim ℓ→∞ S A /ℓ = 0 (see Eq. ( 19)).The combined subsystem A is therefore weakly entangled to the rest of the system, while its two components -one being the mirror image of the other -maintain strong entanglement between them.
The volume-law terms in Eqs. ( 7)-( 10) are evidently generated by the scattering states within the voltage window, with the contribution of each state in Eq. ( 7) being the equivalent of the statistical moment of its corresponding transmission probability.This simple form allows to deduce that the source of the long-range entanglement is the coherence between the reflected part and the transmitted part of each scattered particle, which arrive simultaneously at mirroring sites.In the steady state, the constant particle current renders this strong entanglement a stationary property, and the length scale ℓ mirror determines the amount of entanglement as it is proportional to the number of scattered particles shared by the two subsystems.The voltage bias and the non-trivial scattering constitute necessary and generically-sufficient conditions for the extensive terms in Eqs. ( 7)- (10) to not vanish.
To support our general analytical results, we compared them to numerics for a specific model where the scattering is a result of a single impurity at the site m = 02 .For this model, m 0 = 0 and H scat = ϵ 0 c † 0 c 0 in Eq. ( 4), ϵ 0 being the impurity energy.The scattering matrix for this model yields the transmission probability Good agreement with numerics is manifest in Fig. 2, where, focusing on the aforementioned symmetric case with two intervals of length ℓ, we plotted the scaling with ℓ of all three correlation measures for different ratios of ϵ 0 /η.The analytical results of Eqs. ( 8)-( 10) are plotted for ℓ ≥ 50 (with a constant-in-ℓ additive correction term as the only fitting parameter), as for small values of ℓ there is a considerable contribution from subleading terms beyond the leading volume-law term (an exact analytical result for the most dominant subleading term, which is logarithmic in ℓ, is derived in Ref. [56]).
In Fig. 3 we illustrate a rather counter-intuitive attribute of our results, using the example of the single impurity model.For fixed values of ℓ L and ℓ R , we plot the dependence of the MI and the negativity on the positions of the subsystems, and observe that this dependence is non-monotonic.Indeed, Eqs. ( 8)- (10) state that the long-range correlations are the strongest when the overlap between one subsystem and the mirror image of the other is maximal; if one subsystem is then brought closer to the other, this overlap is reduced and so are the correlations.Fig. 3 again showcases the good agreement of our analytical results with numerics; the apparent slight deviations can be resolved once logarithmic corrections are accounted for [56].

IV. ANALYTICAL METHOD
This section delineates the main steps in the derivation of Eqs. ( 8)- (10), while the discussion of the various technical steps is deferred to the appendices.Subsec.IV A focuses on the derivation of the formulae for the MI and CI (both of which are combinations of subsystem entropies), while Subsec.IV B deals with the derivation of the negativity asymptotics.
The joint starting point of these computations is the two-point correlation function c † j c m for j, m ∈ A. The absence of interactions entails that the states of the total system and its subsystems are Gaussian, and thus entanglement is fully encoded in two-point correlations [43,44,57,58].The correlation function is given explicitly by where u m (k) is the (unnormalized) single-particle wavefunction amplitude at site m of the scattering state associated with momentum k; namely, Eq. ( 12) is derived in Appendix A, where we also discuss how this expression may be simplified when d i ≫ ℓ i .

A. Mutual information and coherent information
The exact asymptotics of the MI and of the CI in Eqs. ( 8) and (10) were obtained through the computation of the Rényi entropies of A L , A R and A. We describe here the main components of this computation, referring the interested reader to Appendix B for full details.
Within the Gaussian steady state, the Rényi entropies of a subsystem X are reduced to functions of C X , the correlation matrix restricted to X ((C X ) jm = c † j c m where j, m ∈ X) [57].Furthermore, these functions admit a simple series expansion, on which our derivation relied.Namely, the Rényi entropies are given by To obtain an analytical expression for the Rényi entropies, it is therefore sufficient to calculate a general expression for moments Tr[(C X ) p ], with p being positive integers.Each such moment can be expressed in the form of a p-dimensional where we defined k 0 = k p .Each sum over m ∈ X in Eq. ( 15) can be rewritten as an integral over a fictitious variable p ] is then expressed as a 2p-dimensional integral.The specific form of this integral depends on the choice of X, and is relatively involved (see Appendix B).We illustrate schematically the way forward by considering the case of the connected subsystems X = A i .For each of these subsystems, Eq. ( 15) can be cast in the general form where k σj = (−1) σj k j , and where the functions p and are independent of ℓ i .The origin of the exponential term in Eq. ( 16) can be traced back to the explicit forms of the wavefunctions u m (k) in Eq. ( 13), which are superpositions of e ±ikm .
An integral of the form of Eq. ( 16) admits a stationary phase approximation in the large-ℓ i limit [29,59,60]; leading-order contributions come only from summands with − → τ = − → σ , with the asymptotics of Eq. ( 16) given by where (−1) . Substituting the specific expressions for the functions f− → τ , − → σ that satisfy Eq. ( 16), we then find that for all positive integers s.A similar treatment was applied to compute the leading term in ∆ℓ i of moments of C A , producing the same result as in Eq. ( 18) up to replacing C Ai with C A and ℓ i with ∆ℓ L + ∆ℓ R .Summing up these contributions in the series expansion of Eq. ( 14) yields the extensive terms for the Rényi entropies of A L , A R and A, The asymptotics in Eq. ( 19) directly lead to Eq. ( 7), while Eqs.( 8) and ( 10) are obtained by taking the limit n → 1 and substituting the resulting von Neumann entropies into the definitions in Eqs. ( 1) and (3).We stress that the universal dependence of the Rényi entropies on the scattering probabilities results from the fact that, at sites m lying outside the scattering region, the wavefunctions u m (k) in Eq. ( 13) are written only in terms of plane waves and scattering amplitudes.

B. Fermionic negativity
Here we outline the principal steps in the derivation of Eq. ( 9), the asymptotic formula for the fermionic negativity; a more detailed account of the computation appears in Appendix C. The derivation relies on the fact that the negativity E can be obtained as the analytic continuation of the Rényi negativities at n = 1, where E n are evaluated at even values of n [43].In analogy to the Rényi entropies, the Rényi negativities E n can be expressed as functions of C A and of a transformed two-point correlation matrix restricted to A [43,44,58].As shown in Appendix C, this expression for E n leads to the following series expansion: Here each C γ is a transformed version of C A , given by where the entries of C A are ordered such that the first ℓ L indices correspond to sites in A L , and the next ℓ R correspond to sites in A R .
Eq. ( 20) reduces the calculation of the Rényi negativities to that of terms of the form Tr C γ1 C γ2 . . .C γp , which, by using Eq. ( 12), may also be written as where the integral is computed over the domain p .The remainder of the calculation is similar in spirit to that of the Rényi entropies: the explicit forms of the wavefunctions from Eq. ( 13) are substituted into Eq.( 22); Eq. ( 22) is rewritten as a 2p-dimensional integral using p fictitious variables; and finally, this 2p-dimensional integral is estimated through a stationary phase approximation (see Appendix C).This process eventually leads to the following result for every positive integer s: Upon summation of the series in Eq. ( 20), we find that the Rényi negativities are given by and the exact asymptotics of the negativity in Eq. ( 9) is obtained once the limit n → 1 is finally taken.

V. DISCUSSION AND OUTLOOK
In this work we derived the exact asymptotics of correlation measures for a nonequilibrium steady state of 1D noninteracting fermions.We have shown that this state hosts extensive long-range entanglement between subsystems that are on opposite sides of a current-conserving scatterer, provided that their distances from it are similar.The volume-law terms of these measures stem from the extensive number of single-particle wavepackets that originate in the high-chemical-potential reservoir, which are split by the scatterer so that they are in a superposition of being found in either one of the mirroring subsystems.The correlation measures thus exhibit a simple and universal dependence on scattering probabilities, allowing to clearly read off the necessary and sufficient conditions for the generation of this strong long-range entanglement.Apart from the requirement that the scatterer be non-trivial, the essential ingredients are the absence of decoherence mechanisms, along with the extensive excess of particles emerging from one of the reservoirs.
We therefore expect the main features of our results to hold in a wide class of nonequilibrium steady states, including those of integrable interacting systems [8,10,61], as well as when the reservoirs are at finite temperatures, and when the scatterer induces particle gain and loss [29].Similar features should also appear in the dynamics following a quench where two decoupled half-infinite chains are prepared with unequal fillings, and the scatterer is suddenly introduced [11,62,63].All of these scenarios offer intriguing prospects for future studies.It would also be interesting to study the interplay of the effects uncovered by this work with decoherence, which could arise due to an integrabilitybreaking impurity or when the system is coupled to Lindblad baths [64][65][66][67], or the interplay of the same effects with the addition of quadratic pairing terms [68], which break charge conservation, to the Hamiltonian of Eq. ( 4).

Appendix A: Two-point correlations
Here we derive Eq. ( 12), which is the general expression for the two-point correlation function c † j c m for sites outside the scattering region, |j| , |m| > m 0 .We consider a long chain with N ≫ 1 sites, where the small scattering region is located at its center; in the end we will take the thermodynamic limit N → ∞.
An annihilation operator c m may be expanded in terms of annihilation operators corresponding to the singleparticle energy eigenstates.We include only extended scattering states in this expansion, neglecting the contribution of localized bound states, since the amplitude of a bound state wavefunction at any site outside the scattering region decays exponentially with the distance of that site from the scatterer.More concretely, we associate an annihilation operator c k,L (c k,R ) with the scattering state of a particle originating in the left (right) reservoir with momentum k > 0 (k < 0).Then, c m may be written as follows: where the wavefunctions u m (k) are given in Eq. ( 13).In the language of scattering state creation operators, the nonequilibrium steady state analyzed in this work is given by with |vac⟩ being the vacuum state.Substituting Eq. (A1) into the definition of the two-point correlation function, we find that, in the thermodynamic limit N → ∞, the correlation function approaches the integral expression of Eq. ( 12).As explained in Sec.IV, the different correlation measures discussed in this work can all be expressed as functions of two-point correlation matrices restricted to the subsystems of interest.That is, we generally consider the terms given by Eq. ( 12) only for sites j, m ∈ A, and correlation measures can scale at most as O(ℓ L + ℓ R ), given the dimensions of the correlation matrices.This, in turn, implies that in the limit d i /ℓ i → ∞ (with d L − d R kept fixed), calculations of correlation measures can be simplified by first neglecting certain terms in the expressions for the matrix elements c † j c m , and only then calculating the appropriate functions of the correlation matrices.In particular, using Eq. ( 13) we observe that in Eq. ( 12) the correlation function is a sum of integrals, where in each integral the integrand is a product of a function of k that is independent of j, m and an exponent of the form exp [iα j,m k], with α j,m ∈ {± (j ± m)}.Then, the Riemann-Lebesgue lemma leads to the conclusion that when |α j,m | ≫ ℓ L , ℓ R , the contribution of the integral is negligible and may be omitted.This entails that when d L , d R ≫ ℓ L , ℓ R we may use the following approximations for the correlation matrix elements: Our analytical results were all derived based on the full expression for c † j c m in Eq. ( 12).As they indicated that the volume-law terms of the different correlation measures depend on d L and d R only through d L − d R , they were compared in Figs.2-3 to numerical results that were computed in the limit d i /ℓ i → ∞, based on the approximated correlation function in Eq. (A3).A comparison to a numerical calculation that relies on the full expression for the correlation function is provided in Appendix D.

Appendix B: Calculation of the Rényi entropies
In this appendix we describe the analytical method used for the computation of the Rényi entropies S n ] for the subsystems X = A L , A R , A. The final results are given in Eq. ( 19).The results for the Rényi entropies lead directly to the asymptotics of the MI and CI (Eqs.( 8) and ( 10)), as explained in Subsec.IV A. In Subsec.IV A we showed that the calculation of Rényi entropies can be reduced to that of the moments Tr[(C X ) p ] for all positive integers p.We now derive the asymptotics of these moments for the subsystems of interest, starting from their integral expression in Eq. ( 15).The analysis is based on the stationary phase approximation (SPA) [59], and is inspired by the analytical methods of Refs.[29,60].

a. Asymptotics of moments for the connected subsystems
We first consider the case X = A R .We begin by introducing the notation W R (x) = x sin x exp 2i m 0 + d R + 1 2 x , and observing that For convenience, we define the notation k aj = (−1) aj k j for a j ∈ {0, 1}, as well as (−1) ⊗p .We use Eqs.( 13) and (B1) to write where Θ(x) is the Heaviside step function, and where we defined

When plugging Eq. (B2) into the expression for Tr[(C A R )
p ] in Eq. ( 15), we will generally get a sum of 2p-dimensional integrals, each of the form ⊗p , and where the function We apply a change of variables and obtain These are the integrals to which we apply the SPA.The SPA allows us to detect which integrals contribute to the leading-order terms of Tr [(C A R ) p ], and to compute their exact contribution to the linear term in ℓ R .From Eq. (B6) it is evident that the answer to the question of whether F( − → τ , − → σ ) has a leading-order contribution is determined by the values of − → τ and − → σ .We now illustrate the method by focusing on two concrete cases for the choice of − → τ and − → σ .
Assuming that τ j = σ j for every j, we find that Applying the SPA to the innermost (2p − 2)-dimensional integral, with respect to the stationary point of the function where we used the fact that the Hessian H at the stationary point yields If, on the other hand, τ j = σ j for every j ≥ 2 but τ 1 ̸ = σ 1 , then Applying the SPA to the innermost integral will again produce a factor proportional to ℓ −p+1 R , yielding only now the remaining phase factor exp [iℓ R k τ1 (ζ 1 + 1)] will eliminate the extensive contribution, such that can have a contribution that is, at most, constant in ℓ R .
From these two examples, it is straightforward to infer the more general rule that an integral Let us now apply this general conclusion to our problem.Substituting Eq. (B2) into Eq.( 15), we obtain Eq. (B8) tells us that the focus on leading-order terms confines F( − → σ , − → σ ) to an integration subdomain where k σj = k σp for all 1 ≤ j ≤ p; this implies that, for the purpose of calculating the leading-order asymptotics of Tr[(C A R ) p ], some of the terms in the full expressions for Ξ aj−1aj k aj−1 , k aj in Eq. (B3) may be a priori discarded, given that for them the k σj−1 = k σj requirement is satisfied only when k σj−1 = k σj = 0. Namely, we may replace and again Ξ 10 (k j−1 , k j ) = Ξ 01 (k j , k j−1 ) * .Note that in the first summand appearing in the expression for Ξ 11 (k j−1 , k j ), the term in the exponent has an opposite sign compared to all other integrals (including Ξ 00 (k j−1 , k j )).Since we have established that expressions of the form F ( − → τ , − → σ ) in Eq. (B4) contribute to the leading order only when − → σ = − → τ , a leading-order contribution to Eq. (B11) will arise from this integral only for − → a = 1 ⊗p , meaning that we can write Applying the SPA as explained above while using the fact that W R (0) = 1, we thus have The derivation for the case

Asymptotics of moments for the disjoint subsystem
We now consider the case X = A. In the summation over sites m in m∈A u m (k j−1 ) u * m (k j ) (the sum that appears in Eq. ( 15)) we will separate mirroring sites from sites which are not mirrored.For concreteness, we assume that where the subsequent generalization is straightforward.We then have We define the function W L (x) = x sin x exp 2i m 0 + d L + 1 2 x .Sums of exponents appearing in Eq. (B16) can be written as integrals: The substitution of Eq. (B16) into the integral expression for Tr[(C A ) p ] in Eq. ( 15) will then yield a sum of integrals of the form where (A j , B j ) ∈ {(∆ℓ L , 0) , (ℓ mirror , ∆ℓ L ) , (∆ℓ R , ℓ L )}.Writing A j = α j ℓ with α j being some fixed ratios, we are interested in the leading-order behavior as ℓ → ∞.Again defining the variables {ζ j } as in Eq. (B5), we arrive at the crucial observation that unless A 1 = A 2 = . . .= A p (and hence also B 1 = B 2 = . . .= B p ), we cannot find in the exponent a (2p − 2)-variable function with a stationary point as before, regardless of the values of − → τ , − → σ .
Leading-order contributions will therefore arise only from terms where A 1 = A 2 = . . .= A p .We can thus conclude that where we defined The first two integrals in Eq. (B19) can be treated in the same way in which the equivalent integrals were treated in the cases of the connected subsystems.What therefore remains to be done is to treat M (p) .We define W(x) = W L (x) e 2ix∆ℓ L in order to simplify the notation.In analogy to Eqs. (B11) and (B12), we may discard terms that have no leading-order contribution to M (p) and write and Ξ 10 (k j−1 , k j ) = Ξ 01 (k j , k j−1 ) * .Applying the SPA through the same procedure as before, while recalling the unitarity of the scattering matrix in Eq. ( 5), we then obtain

Asymptotics of the Rényi entropies
Finally, we use the asymptotics of the moments in Eqs.(B14), (B15) and (B24) to derive the Rényi entropies of the subsystems of interest.In particular, we observe that the terms comprising the series expansion in Eq. ( 14) follow the asymptotic scaling This yields Eq. ( 19), which is true for any d L and d R .
by evaluating E n at even values of n and performing an analytic continuation to n = 1.E n can be written in terms of the restricted correlation matrix C A and a transformed correlation matrix C Ξ , facilitating a significant simplification of the calculation as in the case of the Rényi entropies.We write where the matrices C LR and C RL = (C LR ) † represent two-point correlations between a site in A L and another in A R , and we define where The Rényi negativities can then be written as [43,44,58] We now define the polynomials for any even integer n.Here {z γ } and {z γ } are, respectively, the roots of p n and pn , and they satisfy 2 , so that using the definition of pn we may write the Rényi negativities in Eq. (C5) as Now, if we define the modified correlation matrices one may check that By substituting Eq. (C11) into Eq.(C9), and recognizing that = 1, we arrive at the result which then yields the series expansion of E n reported in Eq. ( 20).As explained in Subsec.IV B, writing this series expansion reduces the calculation to that of terms of the form Tr C γ1 C γ2 . . .C γp , corresponding to the general integral expression in Eq. ( 22).We proceed by applying the SPA to the integrals appearing in Eq. ( 22).For concreteness we again assume that The same argument that led to Eq. (B19) allows us to separate the integral into independent leading-order contributions arising from mirrored and unmirrored sites.Namely, to the linear order in ∆ℓ L , ∆ℓ R and ℓ mirror , we have where we defined (C14) The first two integrals in Eq. (C13) have already been treated using the SPA in Appendix B, as they arise (up to a multiplicative constant) in the calculation of the Rényi entropies (see Eqs. (B14) and (B15)).Assuming for concreteness that k F,L > k F,R , we may therefore write 1 − e 2πiγ j n The asymptotics of both integrals can be estimated using the SPA procedure explained before.When applied to the first integral (which corresponds to an equilibrium scenario), this procedure yields  Using the decomposition of the polynomials p n and pn in Eqs.(C6) and (C7), we arrive at Eq. ( 23), a result that also captures the case k F,L < k F,R (for which an equivalent derivation applies).By summing the series in Eq. ( 20) and taking the limit n → 1, we then obtain the final result for the leading-order asymptotics of the fermionic negativity, given by Eq. ( 9).
-point correlations B. Calculation of the Rényi entropies a. Asymptotics of moments for the connected subsystems b.Asymptotics of moments for the disjoint subsystem c.Asymptotics of the Rényi entropies C. Calculation of the fermionic negativity D. Additional numerical tests References I. INTRODUCTION

Figure 1 .
Figure 1.(a) Schematic sketch of the model: Red circles mark lattice sites in the scattering region, while sites outside this region are marked in blue.Noninteracting reservoirs with different chemical potentials are connected to the two ends of the chain.See the text (Sec.II) for details regarding the notations.Bottom panels: An intuitive picture for the origin of the steady-state entanglement structure.(b) Any incoming wavepacket (black) is split by the scattering region into a transmitted part (dark gray) and a reflected part (light gray), with amplitudes determined by the associated scattering matrix (see Eq. (5)).The transmitted and reflected parts are coherently correlated and thus generate entanglement.(c) Split wavepackets with energies within the voltage window strongly entangle regions that mirror each other with respect to the position of the scattering region.Correlation measures exhibit long-range volume-law scaling, since the number of split wavepackets shared by these mirroring regions is proportional to their length and independent of their spatial separation.

Figure 2 .
Figure 2. The single impurity model: Scaling of (a) the mutual information, (b) the coherent information, and (c) the fermionic negativity between subsystems AL and AR for the symmetric case ℓL = ℓR = ℓ and dL = dR, in the limit di ≫ ℓi.The analytical results of Eqs.(8)-(10) for ℓ ≥ 50 (lines) are compared to numerical results (dots) for different values of the impurity energy ϵ0, with the Fermi momenta fixed at kF,R = π/2 and kF,L = 2π/3.

Figure 3 .
Figure 3.The single impurity model: (a) The mutual information and (b) the fermionic negativity between subsystems AL and AR as a function of their positions relative to the impurity.We fix ℓL = 100 and ℓR = 200, and observe the dependence on dL − dR in the limit di ≫ ℓi.Analytical results (lines) are compared to numerical results (dots) for different values of the impurity energy ϵ0, with the Fermi momenta fixed at kF,R = π/2 and kF,L = 2π/3.Letting ĀL = {m| −m ∈ AL} denote the mirror image of AL, black dashed vertical lines mark the boundaries of the domain where ĀL ⊂ AR, while gray dashed vertical lines mark the boundaries of the domain where ĀL ∩ AR ̸ = ϕ.

k
aj−1 , k aj Θ k aj ∼ ℓ mirror k while for the second integral we findℓ mirror 4π p [k F,R ,k F,L ] p d p k p j=1 Ξ 00 γj (k j−1 , k j ) ∼ ℓ mirror k F,L k F,R

Figure 4 .
Figure 4.The single impurity model: Scaling of (a) the mutual information, (b) the coherent information, and (c) the fermionic negativity between subsystems AL and AR for the symmetric case ℓL = ℓR = ℓ and dL = dR.(d) The mutual information and (e) the fermionic negativity as a function of dL − dR, when fixing ℓL = 100 and ℓR = 200; letting ĀL = {m| − m ∈ AL} denote the mirror image of AL, black dashed vertical lines mark the boundaries of the domain where ĀL ⊂ AR, while gray dashed vertical lines mark the boundaries of the domain where ĀL ∩ AR ̸ = ϕ.In all the panels, results are computed in the limit di ≫ ℓi; analytical results (lines) are compared to numerical results (dots) for different values of the bias ∆k = kF,L − kF,R, with the lower Fermi momentum fixed at kF,R = π/2, and the impurity energy fixed at ϵ0 = η.