Entanglement Measures in a Nonequilibrium Steady State: Exact Results in One Dimension

Entanglement plays a prominent role in the study of condensed matter many-body systems: Entanglement measures not only quantify the possible use of these systems in quantum information protocols, but also shed light on their physics. However, exact analytical results remain scarce, especially for systems out of equilibrium. In this work we examine a paradigmatic one-dimensional fermionic system that consists of a uniform tight-binding chain with an arbitrary scattering region near its center, which is subject to a DC bias voltage at zero temperature. The system is thus held in a current-carrying nonequilibrium steady state, which can nevertheless be described by a pure quantum state. Using a generalization of the Fisher-Hartwig conjecture, we present an exact calculation of the bipartite entanglement entropy of a subsystem with its complement, and show that the scaling of entanglement with the length of the subsystem is highly unusual, containing both a volume-law linear term and a logarithmic term. The linear term is related to imperfect transmission due to scattering, and provides a generalization of the Levitov-Lesovik full counting statistics formula. The logarithmic term arises from the Fermi discontinuities in the distribution function. Our analysis also produces an exact expression for the particle-number-resolved entanglement. We find that although to leading order entanglement equipartition applies, the first term breaking it grows with the size of the subsystem, a novel behavior not observed in previously studied systems. We apply our general results to a concrete model of a tight-binding chain with a single impurity site, and show that the analytical expressions are in good agreement with numerical calculations. The analytical results are further generalized to accommodate the case of multiple scattering regions.

The prospect of promoting our understanding of nonequilibrium quantum systems through the investigation of their entanglement properties is especially intriguing. The study of quantum many-body phenomena out of equilibrium has shown promising progress in recent years, propelled by the development of suitable quantum simulation platforms in cold atom systems [21][22][23][24][25]. But while the theoretical understanding of quantum many-body nonequilibrium has substantially progressed [26,27], rigorous analytical results are still rare. Entanglement measures have established a route for producing such results [28][29][30][31][32][33][34][35][36][37], a route which this work seeks to advance.
Entanglement entropy [1,2] is a measure usually employed to quantify entanglement within a many-body system in a pure state, represented by a density matrix ρ. Given a bipartition of the total system into subsystems A and B, one obtains the reduced density matrix (RDM) of subsystem A by tracing out the degrees of freedom associated with subsystem B, ρ A = Tr B [ρ]. The von-Neumann entanglement entropy (vNEE) is then defined as Additionally, we denote the nth moment of the RDM as and refer to it as the Rényi moment of order n. Note that this is slightly different than the Rényi entropy, S n = 1 1−n ln (Tr [ρ n A ]). The von-Neumann entropy and the Rényi moments are related through Between these two measures of bipartite entanglement, the vNEE constitutes the more rigorous one in and of itself [1,2]. Nevertheless, Rényi entropies may be used to provide lower bounds to the vNEE and to reconstruct the full entanglement spectrum [17,38,39], and are also accessible to direct experimental measurement [40,41]. Of prime importance is the way in which the vNEE scales with the size of the subsystem in question. This scaling law is considered to be a significant classification criterion, distinguishing between typical phases of condensed matter systems [3]. A prominent example is the renowned area law, under which the vNEE scales linearly with the area of the subsystem's boundary: S ∼ cL d−1 , with L being a typical linear dimension of the subsystem, d being the spatial dimension and c being a constant [42]. The area law generally applies to ground states of gapped local Hamiltonians [43,44], as well as to excited eigenstates of many-body-localized systems [14,45,46], and has generated particular interest due to the fact that states obeying the area law admit an efficient tensor-network representation [20,44,[47][48][49]. Ground states of gapless systems (and specifically critical systems) with a finite and sharp Fermi surface tend to violate the area law by a logarithmic correction, S ∼ cL d−1 ln L [6,7,50,51]. On the other hand, a volume-law scaling -i.e., an extensive scaling of the vNEE, S ∼ cL d -is widely observed in highly-excited states of local thermalizing Hamiltonians [52][53][54][55] and in ground states of non-local Hamiltonians [52].
An additional tool in the analysis of many-body entanglement, that has lately come into increased awareness, is the symmetry-or charge-resolved entanglement entropy [56][57][58][59][60]. Given an additive quantity Q = Q A + Q B that is globally conserved in the total system (e.g., particle number), the RDM ρ A derived for a pure eigenstate of the Hamiltonian turns out to be block-diagonal with respect to the eigenbasis of . This suggests that entanglement entropies may be calculated for each block separately, thus resolving the total Rényi moments and vNEE to sums over contributions from symmetry sectors: We note that the definition in Eq. (4) follows Refs. [56,57], while in other works [58][59][60] each symmetry block is normalized by its trace before the resolved moments and entropies are calculated, rendering them measures of entanglement following a projection onto a symmetry sector. These normalized quantities may be straightforwardly derived from their nonnormalized counterparts from Eq. (4) by relying on the fact that Z 1 (Q A ) -which is simply the charge distribution in subsystem A -returns the required trace of each block. Specifically, the post-projection vNEE is given by The resolved quantities in Eq. (4) do not quantify entanglement when used alone, but are more readily calculated, and can be also directly measured in experiments [57,[61][62][63][64]. Symmetry-resolved entanglement represents the internal structure that symmetry imposes on the entanglement spectrum, thus embodying the interplay between conservation laws and entanglement. It has been recently investigated analytically and numerically in various systems, in and out of equilibrium . A common behavior in these systems is entanglement equipartition [58,[67][68][69][70][71][72][73], implying that, to leading order in L, the post-projection vNEE σ (Q A ) is constant across symmetry sectors. The estimation of symmetry-resolved entanglement was shown to yield additional valuable insights, e.g. regarding topological phase transitions [68,78,81] and dissipation in noisy devices [84].
The main result of this work is the exact calculation of an unusual entanglement scaling for the steady state of a one-dimensional fermionic system out of equilibrium, using a generalization of the Fisher-Hartwig conjecture [88]. We study a model of a uniform tight-binding chain containing an arbitrary scattering region at its center, to which a DC bias voltage is applied at zero temperature, thereby leading to a current-carrying steady state. This steady state may be described by a pure eigenstate of the Hamiltonian, with different distributions for scattering states incoming from the left and from the right. This model is relevant to both electronic [89] and cold atom [90] systems.
We report that the subsystem entanglement entropy in this nonequilibrium steady state exhibits a volume-law scaling accompanied by an additive logarithmic correction. More precisely, we find that the vNEE of a subsystem of length L, when located far enough from the scattering region, obeys |t (k)| 2 ln |t (k)| 2 + |r (k)| 2 ln |r (k)| 2    L + C log ln L + C const .
Here |t (k)| 2 and |r (k)| 2 are the transmission and reflection factors (respectively) of the scatterer for a plane wave with momentum k, 0 ≤ k − < k + are the two different Fermi momenta for right-and left-propagating fermions, and C log , C const are constants. The extensive term of the vNEE is thus generated by momentum eigenstates within the bias voltage window, and the contribution of each state is equivalent to the classical mixture entropy of the corresponding transmission probability. The logarithmic term in Eq. (6) is a zero-temperature effect, that arises due to the sharp Fermi-Dirac jumps in the distribution of the plane waves. The coefficient C log is therefore a function of the Fermi momenta, for which we provide an exact expression as well. For C const we obtain an approximate expression, justified under the assumption of a small bias voltage. Furthermore, we expand the result in Eq. (6) by deriving the exact entanglement entropy asymptotics for the case where the chain hosts multiple interspersed scattering regions, finding a similar scaling law. This type of scaling has been previously encountered in rather specific instances, for example in certain excited states, in the ground state of 1D systems with long-range couplings [52,91], and in the diffusive time-averaged state of a 1D interacting system [92,93]. We show how this entanglement scaling can arise generically within a local 1D system in a time-independent state, generalizing previous results in particular cases [29][30][31][32][33].
Our calculation also produces analytical results for the symmetry-resolved entanglement, with respect to the total fermionic charge that is conserved in the system. While the postprojection vNEE exhibits equipartition to leading order as usual, we remarkably find that the first term breaking equipartition may grow with the size of the subsystem, and is antisymmetric in Q A − Q A , the deviation of the charge in A from its mean. To the best of our knowledge, this is the first system found to display such behavior.
The paper is organized as follows: In Sec. 2 we describe the general 1D nonequilibrium model, and derive expressions for the elements of its two-site correlation matrix. In Sec. 3 we review the definition of the generating function that captures all resolved and unresolved moments and entanglement entropies. By employing a generalization of the Fisher-Hartwig conjecture, we analytically derive the exact leading-order asymptotics of the generating function for a large subsystem, along with an approximate subleading correction. We use the result for the generating function to extract quantities relating to charge statistics and the asymptotic scaling of the vNEE in the nonequilibrium steady state. We also address charge-resolved entanglement, and determine the form of the first term breaking entanglement equipartition. In Sec. 4 we apply our general scheme to a specific model, where a single impurity site is responsible for the scattering. We use this concrete example to demonstrate central aspects of our results, and to show that the analytical calculation generally compares favorably with numerics. Sec. 5 details how the analytical results may be generalized to a subsystem on a chain containing multiple but distant scattering regions. In Sec. 6 we discuss our main conclusions and outline possible future directions. Appendix A is dedicated to further technical details of the derivation. Appendix B contains additional plots of the analytical results for the single impurity model, and points out the generic features that apply to the general model.

Model
The general model with which this paper is concerned is that of a long homogeneous onedimensional fermionic tight-binding chain, which contains a finite inhomogeneous region near its middle that induces scattering, as schematically depicted in Fig. 1(a). The single-particle Hamiltonian of such a system can be written as where t > 0 is the hopping amplitude, N + 1 is the number of chain sites (N is assumed to be even), and H scat is the term that will give rise to scattering. H scat should involve only a few chain sites in the vicinity of n = 0, and possibly also hopping terms to additional side-attached sites in this small region. n scat ≥ 0 is an integer such that states |n with |n| > n scat do not appear in the matrix representation of H scat . Subsystem A, for which we will estimate the entanglement measures, includes L contiguous sites (1 L N ) located to one side of the scattering region, and far away from it such that every site n in A obeys |n| − n scat L. For the sake of simplifying the notations we shall focus on the case where the subsystem is to the right of the scattering region, though this is of course an arbitrary choice.
The scattering matrix [94] related to this problem is a function of k, the lattice momentum, S is unitary, and in particular |r L (k)| 2 +|t R (k)| 2 = |r R (k)| 2 +|t L (k)| 2 = 1. For the scattering amplitudes we use the convention that attributes a momentum k > 0 to a state incoming from the left and a momentum −k < 0 to a state incoming from the right, so that S (k) is defined for k > 0. Scattering states constitute the single-particle energy eigenstates for energies |E| < 2t. We assume that the scattering potential does not support a half-bound state, i.e. a non-normalizable solution with energy E = 2t or E = −2t, as is the generic case [95][96][97]. This condition entails that t R (k) , t L (k) → 0 as k → 0, π [98]. The addition of H scat may, however, create bound states in the single-particle energy spectrum, with energies |E| > 2t. If we denote by |ψ (L) k the scattering state related to a wave incoming from the left with momentum k > 0, its form outside the scattering region will be while for a wave incoming from the right with momentum −k < 0, we will have In the many-particle picture, the fermionic creation operator for a site located to the right of the scattering region, n > n scat , can thus be expanded as Figure 1: (a) Schematic illustration of the general lattice model under consideration in Secs. 2-3, with the single-particle Hamiltonian given in Eq. (7). Sites marked in blue belong to the unperturbed parts of the tight-binding chain, while sites marked in red belong to the scattering region, which terminates at the sites n = ±n scat . µ L (µ R ) designates the chemical potential for particles incoming from the left (right). A denotes the subsystem of L contiguous sites with respect to which the calculations of bipartite entanglement in Sec. 3 are performed. (b) Schematic plot of the distributions as functions of lattice momentum, for the scattering states (defined in Eqs. (9)-(10)) and for plane waves (in the region n > n scat , to the right of the scattering region). The case shown is k F,L > k F,R .
where a † k,R creates the scattering state |ψ , and a † k,L creates the scattering state |ψ . If H scat indeed supports bound states, such states would appear as well in the superposition of energy eigenstates that defines |n . We have nevertheless ignored these states in our writing a † n in terms of creation operators of energy eigenstates, due to the localized nature of such bound states, which makes their contribution to the entanglement exponentially small in the distance of subsystem A from the scattering region. Subject to an external constant bias voltage, such a system would arrive at a currentcarrying steady state 1 where the Fermi momentum for waves incoming from the left, k F,L , differs from the Fermi momentum for waves incoming from the right, k F,R . At zero tempera-1 A true steady state is reached only in the case of an infinite chain, N → ∞, which is the limit examined here. One may consider a scenario where the system is prepared using a quench [99], i.e. by connecting at a specific point in time two separate leads with different chemical potentials via the scattering region. Once enough time has passed so that excitations that crossed from one lead to the other have traversed the finite subsystem, the particle fluxes incoming into and outgoing out of the subsystem become balanced, and the entanglement properties of the subsystem relax to time-independent values [36]. ture this steady state is described by a pure many-body state, where |0 is the vacuum state. Taking the limit N → ∞, we replace sums over k with appropriate integrals. The correlation between two sites m, n > n scat for this steady state is then given by where and, in the case where k F,R < k F,L , If instead k F,R > k F,L , we obtain In the case where m, n < −n scat , the result for C mn is similar up to the replacements R ↔ L, τ (k) −→ τ (−k) and h (k) −→ h (−k). As mentioned, subsystem A is assumed to be located to the right of the scattering region, such that n > n scat for every site in A. For further convenience, we denote from now on k − = min {k F,R , k F,L } and k + = max {k F,R , k F,L }. For each momentum k, the distribution n k of the scattering states equals either 0 or 1 by the definition of the steady state. In contrast, the symbol τ (k) may be interpreted as the steady-state distribution of the unperturbed plane waves in the region n > n scat ; a plane wave state within the bias voltage window k − < k < k + will generically be occupied with a fractional distribution factor 0 < n k < 1. This distinction between the distributions of the scattering states and the plane waves is illustrated in Fig. 1(b). As we shall see in Subsec. 3.4, the fractional occupation of the plane waves within the voltage window is the source of the extensive scaling of the subsystem entanglement.
We may conclude that the correlation matrix C is a sum of a Toeplitz matrix, depending on the index difference m − n (the integral in Eq. (13) containing τ (k)), and a Hankel matrix, depending on the sum of indexes, m + n (the integral in Eq. (13) containing h (k)). The Hankel term is negligible for m, n n scat by virtue of the Riemann-Lebesgue lemma, and its decay is generically algebraic, 1 [100]. We can thus write As we elaborate in Sec. 3, all analytical calculations included in this paper rely on the approximation in Eq. (17) of the two-site correlation matrix C. The use of this approximation is justified if the distance between subsystem A and the scattering region is much larger than the length of A, as we further discuss in Subsec. 4.3. There we numerically demonstrate the algebraic decay of the contribution of the Hankel term (neglected within the approximation in Eq. (17)) to the Rényi moments, and also find that this contribution exhibits Friedel oscillations [101,102].

Analytical asymptotics of the entanglement
In this section we present an analytical calculation of the Rényi moments and of the vNEE for subsystem A with respect to its complement. Since in our model the total number of fermions in the lattice is conserved, the Rényi moments and the vNEE may also be resolved with respect to Q A = m∈A a † m a m , the charge in subsystem A. The results we pursue are more conveniently calculated if we start by defining the following generating function [57,58]: where −π < α < π. The calculation of Z n (α) actually encompasses all resolved and unresolved entropies and moments. Rényi moments are given by Z n = Z n (α = 0), from which the vNEE can be extracted through Eq. (3). Furthermore, the charge-resolved Rényi moment is simply a Fourier decomposition of the corresponding generating function [57], In particular, Z 1 (α) is the characteristic function of the charge distribution in A. The chargeresolved vNEE may be subsequently derived through Z n (α) can be written in terms of the eigenvalues {ν l } of 2C A − I L , where I L is the identity matrix of size L, and C A is the two-site correlation matrix (defined as C in Eq. (13)) restricted to subsystem A. More concretely, we may write [57] which enables us to reformulate the calculation as a problem of contour integration in the complex plane [103]: where D L (λ) = det ((λ + 1) I L − 2C A ) and e

Leading asymptotics of the generating function
The immediate consequence of Eq. (17) is that D L (λ) may be approximated as a Toeplitz de- Here we have denoted a definition that may be more compactly packed into the form Importantly, |ν (k)| ≤ 1. The symbol φ (k) as defined in Eq. (23) cannot be cast in the Fisher-Hartwig form [88], contrary to what is required by the well-known (and proven) formulae and theorems of which we are aware [88,104,105] concerning the asymptotics of Toeplitz determinants. A generalized asymptotic formula for the determinant of a Toeplitz matrix generated by a piecewisecontinuous symbol was conjectured in Refs. [106,107]. According to this formula, for L 1 we have where the ellipses represent terms of lower order in L.
Plugging the asymptotic form in Eq. (26) into the integral expression in Eq. (22), we obtain where we have defined The term I lin (n, α) is derived in a straightforward manner, while a detailed derivation of the term I log (n, α) appears in Appendix A.1.
The linear term in L appearing in Eq. (27) counts the filled momentum (plane wave) states. The states with −k F,R < k < k − are all filled with probability 1, while if k F,R < k F,L , the states with k − < k < k + are filled with probability |t L (k)| 2 and empty with probability 1 − |t L (k)| 2 (or vice versa in the case where k F,L < k F,R ). This distribution is schematically presented in Fig. 1(b). Since any interval δk includes (L/2π) δk states, exp (I lin (n, α) L) is simply the product of the moments that arise from the individual filled states. The term exp (I lin (n, α) L) can thus be interpreted as a generalization of the generating function for the full counting statistics in the case of a transmission factor that is constant in k, cf. Eq. (20) in Ref. [99] (the Levitov-Lesovik formula).
The logarithmic term in L appearing in Eq. (27) is a result of the Fermi discontinuities at k = k ± and k = −k F,R , featured in Fig. 1(b). In particular, since Q n (1, α) = 0 by definition, the contribution from k = k ± vanishes when ν (k ± ) → ∓1 (respectively), in accordance with the disappearance of the respective jump discontinuity of the symbol φ (k). That is also the case when the bias voltage is larger than the bandwidth, such that k F,L = 0 and k F,R = π or vice versa: since |r R (k)| → 1 as k → 0, π, the symbol φ (k) is continuous at any value of k, and thus the logarithmic term vanishes from Eq. (27).

Subleading asymptotics of the generating function
In order to incorporate further subleading terms in the analytical asymptotics of ln Z n (α), we now approximate the symbol φ (k) in Eq. (23) to be piecewise-constant, such that it will fit the Fisher-Hartwig form [88]. The underlying assumption is that under a small bias such . This requires the transmission and reflection factors to change slowly near k = k 0 . For this purpose we define the following approximate symbol:φ where we denoted ν 0 = ν (k 0 ).
Let us denote byD L (λ) the determinant of the Toeplitz matrix generated by the symbol in Eq. (29). Using the Fisher-Hartwig formula [88,103] for the asymptotics ofD L (λ) and substituting it into the integral expression (22) for the generating function will yield for L 1 an expression of the form Note that while the asymptotics of ln Z n (α) in Eq. (27) was estimated up to a linear term and a logarithmic term in L without an approximation of the Toeplitz symbol, here the approximationφ (k) yields an additional term which is independent of L. Compared to the exact expressions for the linear and logarithmic terms, the error of the terms obtained from the approximate symbol scale as as ∆k → 0.
Relying on the Fisher-Hartwig formula, we obtain that where we have defined We detail the calculation ofĨ const (n, α) in Appendix A.2. In total, the approximate asymptotic expression for the generating function is Note that, generically, the terms in I log (n, α) andĨ const (n, α) which stem from the partial transmission effects for k − < k < k + do not vanish when we fix L and simply take the limit ∆k → 0; we must first take ν (k) → ±1 in the interval k − < k < k + in order for them to vanish. This, however, is in compliance with the fact that the problem is examined at the limit of large L, which thus constitutes the largest length scale of the problem (other than the length scales that are assumed infinite, i.e. the length of the full chain and the distance of A from the scattering region). The limits L → ∞ and ∆k → 0 do not commute; in other words, our expressions assume ∆k 1/L, and hence taking ∆k → 0 in them does not eliminate the contributions coming from the partial transmission within the voltage window ∆k.
A more problematic feature of the analytical result is that for any n, the function Q n (ν, α) (which appears in the expressions for both I log (n, α) andĨ const (n, α)) is singular at ν = 0 and α = ±π. Moreover, it may be shown that This deems the analytical expression for Z n (α) to be a non-integrable function of α over [−π, π] if either ν (k + ) = 0 or ν (k − ) = 0. Numerical results do not exhibit this kind of divergence, and so this property of the analytical result does not capture the true behavior of the generating function. The singularity of the function e (α) n (1, ν) (defined right after Eq. (22)) at ν = 0, α = ±π is what brings about this difficulty, as small shifts of ν (k ± ) become crucial when either one is near ν = 0.
We note that our previous work [68], which had discussed a situation where the Toeplitz symbol may indeed be cast in the Fisher-Hartwig form, established that corrections to the approximation of ln Z n (α) using the Fisher-Hartwig formula decay less rapidly with L as α nears ±π. At α = ±π these corrections eventually become as important as the terms in Eq. (30), causing a considerable deviation from exact numerical results if one does not include the corrections [67,68,70]. Although there is no known expression for these corrections when the Toeplitz symbol does not fit the Fisher-Hartwig form, we expect them to eliminate the divergence at α = ±π observed in this case. While the divergence of the generating function prevents us from obtaining analytical results for charge-resolved quantities in cases where either ν (k + ) or ν (k − ) exactly vanish, we have found that it has little effect whenever ν (k ± ) are finite, even when they are small, as demonstrated in Subsec. 4.2.
We may also estimate the deviation of the generating function for the nonequilibrium steady state from that of the ground state in the equilibrium case. In Refs. [67,68] it was shown that the leading-order asymptotics of the generating function for the ground state of a homogeneous tight-binding chain filled up to k = ±k 0 is given by 2 Subtracting this from the nonequilibrium result, we obtain

Charge statistics
An important special case of the symmetry-resolved Rényi moments is that of n = 1, since Z 1 (Q A ) constitutes the charge distribution in subsystem A, and an expansion of ln Z 1 (α) in powers of α gives its moments. Indeed, we may write where is the shift in the mean charge, and is the shift in the charge variance. Here γ E ≈ 0.577 is the Euler-Mascheroni constant [108]. A derivation of Eqs. (39) and (40) Analogously, we may define a generalized quantity Q A n = −i∂ α ln Z n (α) | α=0 , designating the mean of the "charge distribution" whose characteristic function is Z n (α). By taking the derivative of Eq. (34), this generalized mean charge is found to be In similar fashion, one can obtain an analytical expression for the corresponding variance by calculating (∆Q A ) 2 n = −∂ 2 α ln Z n (α) | α=0 , and in particular find that generically it scales linearly with L.

Unresolved entanglement
By setting α = 0 in the generating function from Eq. (34), we obtain analytical expressions for the unresolved entanglement measures we wished to estimate. Rényi moments and entropies are directly accessible in this manner, while the vNEE is extracted through its relation to the Rényi moments, per Eq. (3). In order to conveniently present the resultant asymptotics for the vNEE, we define the following functions: where we have introduced the constant − e −z 12z dz ≈ 0.1399. We then have for the vNEE the following result: where C lin and C log are both exact and are given by and and C const is the approximate constant correction (valid, as before, when ∆k is small enough), This is derived with further details in Appendix A.4. Eq. (44) was already highlighted in Sec. 1 (where it is featured as Eq. (6)) as a central result of this work.
The form of C lin in Eq. (45) is especially illuminating. The extensive term of the entanglement entropy arises from the integration of a classical mixture entropy with respect to the reflection and transmission probabilities. It also highlights the two crucial ingredients that produce the linear leading term: the nonequilibrium setting brought about by the bias voltage, which ensures that k + = k − and therefore that the integral does not trivially vanish; and a scattering potential that generates imperfect transmission, seeing that a unity transmission factor will cause the integrand in Eq. (45) to vanish for all k. Both conditions must apply in order for the leading term of S to be extensive (for related treatments of particular time-dependent impurity setups, see Refs. [29,32]).

Charge-resolved entanglement
An exact computation of the charge-resolved Rényi moments based on the analytical asymptotics of the generating function in Eq. (34) requires carrying out the integration in Eq. (19), which cannot itself be performed analytically. For L 1, a useful approximation is obtained by expanding ln Z n (α) in powers of α up to second order and replacing the integration limits in Eq. (19) by ±∞. This leads to an approximate Gaussian form of the charge-resolved nth Rényi moment: The above approximation holds since for large L, (∆Q A ) 2 n scales linearly with L (as mentioned in Subsec. 3.3), and consequently Z n (α) decays rapidly away from α = 0 (this is analogous to the central limit theorem).
The Gaussian approximation allows us to analytically examine the question of entanglement equipartition. By plugging Eq. (48) into Eq. (5), we obtain the following expression for the vNEE after a projective charge measurement: Relying on Eqs. (40) and (42), we note that both (∆Q A ) 2 and ∂ n Q A n | n=1 scale linearly with L to leading order. In particular, the approximation in Eq. (49) is expected to be valid for values of Eq. (49) entails that, to leading (linear in L) order, entanglement is spread equally among charge sectors with σ (Q A ) ≈ S, but also that this equipartition may be broken by a term up to order O √ L . To the best of our knowledge, this is the first calculation of symmetry-resolved entanglement entropy showing a term breaking equipartition that grows with the size of the subsystem in question [67][68][69][70][71][72][73].
The fact that the first term breaking entanglement equipartition is odd with respect to Q A − Q A is noteworthy as well. Previous works that explicitly calculated the first equipartition-breaking term in different equilibrium and nonequilibrium models have always found it to be an even function of the deviation from the mean charge [67][68][69][70][71][72][73]. For L 1 this odd term is given by which suggests that, for a small nonzero bias voltage, the post-measurement vNEEs of two charge sectors Q A = q 1 and Q A = q 2 approximately differ by In the cases of either zero bias voltage (∆k = 0) or perfect transmission or reflection (t L (k) = 1 or t L (k) = 0, respectively, for all k), we have ∂ n Q A n = 0, so the equipartitionbreaking term that displays both these novel features vanishes. We stress that these features are unique even with respect to the case studied in Ref. [73], where entanglement equipartition is examined for a nonequilibrium steady state created following a global quench. The main distinction between the steady state in Ref. [73] and the steady state we investigated here is that the latter is a state with a net current that is partially transmitted by the scatterer.

The single impurity model
We consider a concrete example of the general model discussed above, by setting an on-site energy cost for the middle site of the chain. The single-particle Hamiltonian in Eq. (7) becomes where 0 ∈ R. The scattering states which constitute solutions for the single-particle problem provide the following transmission and reflection coefficients: There is also a bound state with energy E > 2t (E < −2t) for 0 > 0 ( 0 < 0); however, as noted above, its contribution is negligible in the limit considered. The generating function Z n (α) thus depends on the parameters k ± and 0 /t, along with the more explicit dependence on α and n. Appendix B illustrates how the coefficients that define the analytical asymptotic expression for ln Z n (α) in Eq. (34) vary with these parameters. In this section we first focus on measures of entanglement extracted from the generating function for the single impurity model, with a comparison of our analytical results to numerics (Subsecs. 4.1 and 4.2). We then use numerics for this model to discuss more generally the accuracy of the calculation of the generating function (Subsec. 4.3).

Unresolved von-Neumann entanglement entropy
As already stressed, the entanglement between subsystem A and its complement is most rigorously quantified by the vNEE. In Fig. 3 we plot the dependence on the model parameters of the asymptotic scaling coefficients of the vNEE from Eq. (44). It can be seen that the coefficients are nonmonotonic in 0 /t. Notably, the integral form of C lin in Eq. (45) suggests that the largest contribution to the leading extensive term of S comes form momentum states where |t L (k)| 2 ≈ |r R (k)| 2 ; this, in turn, implies that C lin should peak at a value of 0 /t such that |t L (k 0 )| 2 ≈ |r R (k 0 )| 2 , as is evident in Fig. 3(a).
Another noteworthy detail is that Fig. 3(b) exemplifies the non-continuous nature of the asymptotic scaling coefficients that depend on the Fermi discontinuities. Indeed, in a homogeneous chain ( 0 = 0) we have C log = 1/3 for any 0 < k 0 < π, but Fig. 3(b) shows that by fixing k F,L = π first and taking the limit 0 /t → 0 later, we obtain C log → 1/6. This is because for any 0 = 0 there is no Fermi discontinuity at k F,L = π, but one is created at 0 = 0. One must therefore be careful when taking limits that either create or destroy jumps in the distribution.
To corroborate these analytical results, the vNEE was also extracted through Eq. (3) from a numerical calculation of the Rényi moments Z n . The latter is performed using the exact expression for the generating function in Eq. (21), where the restricted correlation matrix C A is approximated according to Eq. (17), so that effects of a finite distance between subsystem A and the impurity site are neglected. Figs. 3(d)-(e) feature a comparison of the numerical result for the vNEE with our analytical calculation, confirming good agreement between them. As attested by Fig. 3(d), this is true even for a bias voltage that is not relatively small, such that ∆k ≈ 1. Recall that the assumption ∆k 1 was required only for justifying the approximation of C const , while both leading terms of the asymptotics are exact regardless of it.

Charge-resolved entanglement
Next, we studied the symmetry resolution of the Rényi moments Z 1 and Z 2 and of the vNEE. This was done by extracting the symmetry-resolved quantities from both the analytical calculation (Eq. (34)) and the numerical calculation of the generating function Z n (α), relying on Eqs. (19) and (20). The numerical estimation of Z n (α) was obtained using Eq. (21), again using the approximation in Eq. (17) for the correlation matrix, which assumes an infinite distance between subsystem A and the impurity.
The results are presented in Fig. 4, where it is evident that the numerical results (naturally sampled at integer values of Q A ) accurately fit the analytical results near the mean charges of the distributions, given by Q A n from Eq. (42) (with n = 1 for S (Q A ) and Z 1 (Q A ), further simplified in Eq. (39), and with n = 2 for Z 2 (Q A )). The plots in Fig. 4 are all centered around the integer charge that is the nearest to the corresponding mean charge, given by where is the floor function ( x is the nearest integer to x from below) and is the ceiling function ( x is the nearest integer to x from above). The charge-resolved quantities indeed reach their maximal value near the mean charge Q A n , but due to their slight deviation from a Gaussian form (since L is large but finite) the charge sector is not necessarily the sector where they peak. Note that the mean charge Q A n varies with the model parameters 0 /t and k ± .
A conspicuous property of the resolved quantities is that their distribution among charge sectors becomes wider as the model parameters approach values such that |t L (k)| 2 ≈ |r R (k)| 2 for k − < k < k + . This is manifested in the analytical results most simply for Z 1 (Q A ), since  Fig. 4 we fixed k − = π/2 and ∆k = 0.1, so this condition is equivalent there to 0 ≈ 2t. The exact point 0 = 2t cannot be examined analytically due to a non-integrable divergence of the generating function, as explained in Subsec. 3.2. Nevertheless, the results for 0 = 1.99t in Fig. 4 indicate that, even at points very close to the specific point where the divergence occurs, this divergence does not cause any discernible deviation of the analytical calculation from numerical results.
Additionally, we examined the post-projection charge-resolved vNEE, σ (Q A ), for the single impurity model. Since the analytical result of Eq. (49) relies on the Gaussian approximation of the generating function, it was natural to test whether it properly captured the behavior of the charge-resolved measures that were extracted from the more accurate analytical form of the generating function, given by Eq. (34). In Fig. 5 we present a comparison between the deviation from entanglement equipartition of the full analytical result and the leading-order term estimated in Eq. (50), which is linear in Q A − Q A . From Fig. 5(a) it is evident that the linear breaking of entanglement equipartition holds up to |Q A − Q A | ≈ ∆Q A , confirming that the equipartition-breaking term scales as √ L for large L. Fig. 5(b) affirms that for large L, the slope of this linear term becomes independent of L. Fig. 5(b) also includes fully numerical estimations (using the numerical calculation of Z n (α) and Eqs. (5), (19) and (20)) of the change in σ (Q A ) between adjacent (integer valued) charge sectors near Q A , for reasonable subsystem sizes. These numerical results nicely follow the trend of their analytical counterparts, once again attesting to the validity of the latter.

Accuracy of the generating function calculation
With the generating function Z n (α) being the basis for all the analytical calculations discussed in this paper, a test of the accuracy of its calculation across the parameter space is required. In Fig. 6 the analytical estimation of ln Z n (α) (denoted as ln Z (ana) n (α)) for the single impurity model is compared to a numerical calculation of ln Z n (α) (denoted as ln Z (num) n (α)). Here the numerical calculation once again neglects the effects of a finite distance between subsystem A and the impurity (whose effects will be examined later on), relying on Eqs. (17) and (21). The comparison indicates good agreement between analytical and numerical results for values of α far enough from α = ±π. As explained above, the fact that this agreement breaks down as α approaches ±π is a well-known trait of the leading-order approximation that stems from the Fisher-Hartwig conjecture [67,68,70].
Importantly, Fig. 6 illustrates that for fixed values of k ± , divergences in the absolute deviation between analytical and numerical results may occur at α = ±π, though they are typically rare within the parameter space. Such singularities appear if either ν (k − ) = 0, ν (k + ) = 0 or ν 0 = 0, where ln Z Eq. (35)), or if one of the numerically calculated eigenvalues ν l vanishes, thus setting the numerical result to Z (num) n (α = ±π) = 0 according to Eq. (21). In Fig. 6(b), for example, singularities at α = π may be detected for values of 0 near 2t (divergence of ln Z (ana) n (α)) and near 1.5t, 2.5t (points where Z (num) n (α) vanishes; the locations of these points in the parameter space varies with L). However, when varying k ± and fixing all other parameters, ln Z (num) n (α) oscillates rapidly as a function of ∆k, with periodicity 2π/L. These oscillations are not reflected in the analytical result, causing its deviation from ln Z (num) n (α) to oscillate rapidly as well, as seen in Fig. 6(c). Fig. 6(c) more specifically indicates that for α = ±π, there exists a regime of values of k ± where in each period of the oscillation there is a singularity of the deviation. Each such singularity corresponds to a vanishing numerically calculated eigenvalue ν l , i.e. to a zero of Z (num) n (α = ±π). Finally, we address the effect of a finite distance between subsystem A and the impurity -i.e., of the Hankel term in Eq. (13) that was omitted from Eq. (17), and was heretofore disregarded. Although we were not able to incorporate the effect of the Hankel term into our analytical calculation, the concrete example of the single impurity model allows us to examine numerically the dependence of the results on d, the distance of A from the impurity at the origin. More precisely, d is taken to be the location of the leftmost site in subsystem A, such that A includes the sites n = d, . . . , d + L − 1. Fig. 7 shows the comparison between two numerical calculations of the generating function: the one extracted from the approximate form of the correlation matrix in Eq. (17) (denoted as Z (∞) n (α)), which was used in the comparison to the analytical results, and the one that relies on the full form of the correlation matrix in Eq. (13) (denoted as Z (d) n (α)). One may observe that the additional term that accounts for finite d effects on ln Z n (α) oscillates as a function of d, a manifestation of Friedel oscillations [101,102]. The typical wavenumber of these oscillations is 2k F,R , and is therefore independent of k F,L , as should be expected from the fact that the d-dependent term in Eq. (13) is independent of k F,L as well (this holds for a subsystem to the right of the scattering region; for a subsystem on the left, the roles of k F,L and k F,R are switched). The numerical results also verify that the deviation of ln Z  [101]. This observation dovetails with the aforementioned projection of an algebraic decay of the Hankel term in the correlation matrix.
The dependence of the generating function on d may be seen as an effect of boundary conditions, referring here to the boundary of subsystem A. Finite d effects are therefore expected not to be reflected in the leading, linear in L term of ln Z n (α), since this term represents an extensive property of A. We have verified numerically that the Hankel contribution indeed has no extensive effect. The logarithmic term of ln Z n (α), in contrast, does generally depend on the boundary conditions of A (cf. Refs. [51,109]), and so for finite d the Hankel contribution may be proportional to ln L, though such a logarithmic dependence is in practice hard to ascertain numerically for accessible subsystem sizes. Regardless of this, considering that in the previous calculations we kept terms which are constant in L, and that the Hankel contribution is proportional to d −1 provided that d L, neglecting the d dependence is certainly justified in the regime d L.

Generalization to multiple scatterers
In the following section we rely on the analytical results for the model described in Sec. 2 in order to derive corresponding results for a more general scenario, where the tight-binding chain contains several different scattering regions rather than just a single one. The necessary foundation is the description of the combined scattering effects of two scatterers, with a distance of sites between them. We assume throughout this section that is considerably larger than the Fermi wavelengths, (k F,L ) −1 and (k F,R ) −1 . We mark the left scatterer with I, and the right one with II. Let us associate a unitary scattering matrix with each scattering region, where i = I, II. If subsystem A is situated to the same side of both scattering regions, as depicted in Fig. 8(a), we may treat them as a single scatterer with appropriate reflection and transmission amplitudes [89]. For a wave incoming from the left, these are given by and for a wave incoming from the right the amplitudes r (I+II) R and t (I+II) R are given by the same expressions, up to exchanging R ↔ L, I ↔ II and multiplying by an overall phase (the notation emphasizing the dependence of the scattering amplitudes on k has been omitted for brevity).
As described in Sec. 3, the correlation matrix C A is the basis for our calculation of entanglement measures. By the same line of argument of Sec. 2, it can be written as a sum of a Toeplitz matrix and a Hankel matrix, where the Hankel matrix can be neglected assuming A is far enough from the combined scatterer, as before. For the estimation of the Toeplitz term, , denotes the number of sites between scattering regions I and II, which is assumed to be much larger than the Fermi wavelengths.
we define the following incoherent scattering probabilities [89]: The probabilities |r L (k)| 2 and |t L (k)| 2 differ from the expressions r , respectively, under integrals over k. This finally allows us to write the two-point correlation matrix of A according to Eq. (17), where the incoherent scattering probabilities stand in for the original scattering probabilities of a single scattering region. This scheme can be readily extended to treat multiple scattering regions that are all located on the same side of subsystem A (as illustrated in Fig. 8(b)) as a single combined scattering region, and therefore one may use the analytical results of Sec. 3 to estimate the various entanglement measures discussed there.
Considering the case where scattering regions are located on both sides of subsystem A (depicted generally in Fig. 8(d)), this scheme of calculating combined scattering probabilities has reduced the problem to that of two scattering regions -region I to the left of A, and region II to its right, as illustrated in Fig. 8(c). The distances from the edges of A to the scatterers are assumed to be much larger than L, the number of sites in A. Using the assumption of the distance between scatterers I and II being much larger than the Fermi wavelengths, we arrive at the following approximation of the correlation matrix in that region: Here we have defined where and ν II (k) is defined similarly, up to replacing t (I) The derivation of Eq. (58) is detailed in Appendix A.5.
The analytical method of Subsec. 3.1 can now be applied to obtain the exact form of the two leading terms in the asymptotic expression for the generating function defined in Eq. (18): An approximation for the first subleading correction to this asymptotics, which is independent of L, may be derived by following the analytical method of Subsec. 3.2. The various entanglement measures that were discussed in the context of a single scattering region can now be extracted from Eq. (61). We report in particular the result for the unresolved vNEE in subsystem A between the two scatterers. For convenience, we introduce the notations The vNEE is given by the asymptotic form where the function q (p) was defined in Eq. (43). Note that if we require transmission and reflection factors to be constant functions of k, we recreate Eq. (26) of Ref. [31]. We have therefore generalized the scenario discussed in Ref. [31], where the subsystem lies on a tight binding chain coupled to two macroscopic leads with k-independent hybridization factors.

Conclusions and outlook
While the exact physical description of nonequilibrium many-body states remains a coveted yet elusive goal, entanglement measures continue to facilitate incremental progress toward its achievement. In this work we sought to exactly quantify the steady-state entanglement of a paradigmatic 1D lattice model, namely a homogeneous tight-binding chain interrupted by an arbitrary number of scattering regions, held under a bias voltage at zero temperature. For this purpose we employed the generalized Fisher-Hartwig conjecture to calculate bipartite entanglement measures for a subsystem located far away from the scatterers. A central result of our work is given in Eq. (44). The von-Neumann entanglement entropy was shown to scale extensively with the size of the subsystem in question, with an additive logarithmic correction that arises from the sharp jumps in the energy distribution. While in the ground state such a scaling law is considered exotic and requires long range couplings [52,91], the class of steady states we investigated, which are excited eigenstates of the Hamiltonian with a Fermi-discontinuous distribution of the excitations [52], exhibits it generically. This suggests that a scaling law of the form of Eq. (44) should be much more common in nonequilibrium setups, rendering it a strong signature of the unique properties that distinguish steady states in and out of equilibrium. More precisely, this entanglement scaling law should in general be observed in steady states where the single-particle energy distribution features partially occupied states (of non-vanishing measure) and Fermi discontinuities. Refs. [30,31] have indeed found such scaling of the vNEE for nonequilibrium steady states in some particular impurity setups. We expect this behavior to apply also to other currentcarrying impurity models at zero temperature, including systems that contain massless Dirac fermions 4 or host dissipative defects [110,111].
Notably, the form of the extensive term of the vNEE encapsulates the basic elements defining the steady state. It arises from scattering states within the energy window between the two different chemical potentials, as scattered particles traverse the subsystem from side to side, and thereby entangle its entire bulk to the rest of the chain. The assumption of a large subsystem allows us to disregard boundary effects and attribute classical probabilities to the scattering processes, and therefore the contribution of each mode to the entanglement is equivalent to a classical mixture entropy. As a consequence, the extensive term of the vNEE vanishes either in the absence of a bias voltage (i.e., in equilibrium), or if the scattering region is trivial. In this sense, the model studied here can be seen as a minimal model for producing such scaling of the steady state entanglement. The picture of entanglement as a result of partial occupation of momentum states due to scattering is also reflected in the linear term of ln Z n (α) in Eq. (27), which generalizes the known full counting statistics formula found by Levitov and Lesovik [99,112,113].
The exact expression for the vNEE is only one out of the comprehensive set of results presented in this paper, all encoded in the asymptotics of the generating function Z n (α) in Eq. (34). These results include Rényi moments (from which Rényi entropies are readily obtained), statistical charge properties, and charge-resolved moments and entanglement measures. We particularly emphasize the novelty of the result in Eq. (49) for the vNEE following a projective charge measurement σ (Q A ). It implies that to leading order, entanglement is equally distributed across charge sectors, as has been established for the great majority of models for which symmetry-resolved entanglement was studied. However, the term breaking entanglement equipartition grows with the subsystem size L, as σ (Q A ) − S = O √ L , for charge sectors within the standard deviation from the mean charge. The leading equipartitionbreaking term was also found to be anti-symmetric in To the best of our knowledge, the model studied in this paper is the first where σ (Q A ) exhibits these properties. A natural question for future research is therefore whether this unique behavior can be exclusively ascribed to the partially-transmitted current carried by the nonequilibrium steady state, as it was not witnessed in a different nonequilibrium model where the net current had been absent [73]. We additionally note that we expect the new behaviors uncovered in this work to hold in the presence of interactions, although further investigation is required to establish this claim.
Our analytical results were tested against numerics in a model where a single impurity site serves as the scatterer, and were shown to compare to them favorably. We additionally used numerics to observe effects of a finite distance between the subsystem and the scattering region, detecting signatures of Friedel oscillations and confirming that the effects are indeed negligible when that distance is large enough. We note that by using various proposed protocols for the measurement of (resolved and unresolved) entanglement measures [61][62][63][64]84], our results may be experimentally tested in setups based on cold atom or electronic systems [89,90].
The road toward a deeper understanding of nonequilibrium many-body physics is still riddled with unanswered questions. Hopefully the exact results presented in this paper can serve as building blocks for the future description of richer and more intricate phenomena out of equilibrium, such as effects of interactions, disorder, localization or external driving [53,114,115], entanglement phase transitions [92,[116][117][118], and transport through mesoscopic systems [119,120].
where we have employed integration by parts. Using the property from Eq. (A.1) and taking the limit ε, δ → 0 + , we obtain and therefore Here we invoked the notation from Eq. (28), and the term written explicitly was obtained by carrying out the integration through the change of variables u = ln x−1 x+1 .

A.2 Subleading term in ln Z n (α)
We detail the calculation of the expression in Eq. (32) for the termĨ const (n, α), which is the approximate subleading term in the asymptotics of ln Z n (α) in Eq. (34). The approximate piecewise-constant symbolφ (k) of Eq. (29) can be written in a Fisher-Hartwig form [88]: , and where we also defined and According to the Fisher-Hartwig conjecture, the asymptotics of the Toeplitz determinant D L (λ), arising from this symbol, is given by [88,103] where G (x) is the Barnes G-function [121], which obeys in particular The required subleading termĨ const (n, α) in Eq. (30) is therefore given by the integral Let us now denote ln e ik j − e ikm 2β j βm d dλ e (α) n (1 + ε, λ) dλ.

A.3 Expansion of the generating function for n = 1
Here we explicitly calculate the terms of order α and α 2 in the power series expansion of ln Z 1 (α). We first note that e (α) 18) and therefore the linear term in Eq. (27) obeys For the logarithmic and constant terms we use the fact that Q 1 (ν, α) may be calculated explicitly for any −1 ≤ ν ≤ 1. Indeed, the change of variables u = x−ν 1−x allows us to write and then solve the integral using complex contour integration of the function ln 2 z/ (1 + z) cos α 2 + i (z + ν) sin α 2 (1 + z). We eventually obtain where the logarithm should be interpreted as belonging to the principal branch, |Im ln z| < π. Expanding in powers of α, we arrive at Finally, we estimate the contributions of terms of the form Υ 1 (ν, α) that appear in the expression for the constant term in Eq. (32). Up to order α 2 we have, according to Eqs. (33) and (A.18), Employing a change of variables ζ = ln x−ν 1−x and the useful formula [103] ln one ends up with (A.25) Switching the order of integration, the O (α) term vanishes trivially due to the integrand being odd with respect to ζ. Carrying out the integration of the rest, we arrive at t(e t −1) dt is the Euler-Mascheroni constant [108]. Adding up the different terms, we conclude that the expansion up to order α 2 of the nonequilibrium deviation of the generating function (Eq. (37)) for n = 1 is given by Applying Eq. (25), the expansion may be also expressed as

A.4 The von-Neumman entanglement entropy
We present here details of the derivation of the asymptotic form for the unresolved vNEE in Eq. (44). From Eqs. (3) and (34), along with the fact that Z 1 = 1 by definition, we draw the following relation: The explicit expression for C lin in Eq. (45) is then straightforward to obtain. Calculating analytically the derivatives of ∂ n I log (n, 0) and ∂ nĨconst (n, 0), however, requires a more subtle analysis. The functions Q n (ν, 0) and Υ n (ν, 0) that appear in those terms (defined in Eqs. (28) and (33)) have integral definitions with integrands that depend on n, and the integrals must be rewritten before one can estimate the required derivatives simply by taking the derivatives of the integrands. Both Q n (ν, 0) and Υ n (ν, 0) are defined as integrals over the interval [ν, 1]. By splitting it into the intervals ν, 1+ν 2 and 1+ν 2 , 1 , and by changing variables to u = x−ν 1−x within the former and u = 1−x x−ν within the latter, we arrive at the following expressions:

A.5 Correlation matrix in a subsystem between two scatterers
We derive the approximate form appearing in Eqs. (58) and (59) for the correlation matrix of a subsystem situated between two scattering regions -region I on its left and region II on its right -assuming that the edges of the subsystem are far away from both regions. Energy eigenstates are given by scattering states, which for a wave incoming from the left take the form which is justified given that the difference between the two expressions oscillates as a function of k with a frequency of 2 , and thus after integration its contribution decays for that is large with respect to the Fermi wavelengths, (k F,L ) −1 and (k F,R ) −1 . We therefore obtain This result holds true regardless of the direction of the bias voltage. Once we explicitly distinguish between the two possible directions, we finally obtain the expression in Eqs. (58) and (59).

B Additional plots for the single impurity model
Focusing on the single impurity model defined in Sec. 4, the plots presented here illustrate the parameter dependence of the coefficients in the analytical asymptotics of the generating function (Eq. (34)). Fig. 9 demonstrates how, as n increases, the dependence of Re ln Z n (α) on α becomes flatter, which amounts to a narrower distribution of the corresponding chargeresolved Rényi moment about its peak. Figs. 10-11 depict the dependence of the coefficients on  27), which is independent of n. and 11(b),(e) the divergences of I log for α = π are related to points where either ν (k − ) = 0 or ν (k + ) = 0. In Figs. 10(c),(f) the coefficientĨ const is seen to diverge for α = π and ν 0 = 0. Note that in Fig. 11, the values of I log andĨ const at the limit ∆k → 0 are fictitious since, as explained in Subsec. 3.2, throughout our calculations we implicitly assume ∆k 1/L.