Entanglement and thermokinetic uncertainty relations in coherent mesoscopic transport

A deeper understanding of the differences between quantum and classical dynamics promises great potential for emerging technologies. Nevertheless, some aspects remain poorly understood, particularly concerning the role of quantum coherence in open quantum systems. On the one hand, coherence leads to entanglement and even nonlocality. On the other, it may lead to a suppression of fluctuations, causing violations of thermo-kinetic uncertainty relations (TUR and KUR) that are valid for classical processes. These represent two different manifestations of coherence, one depending only on the state of the system (static) and one depending on two-time correlation functions (dynamical). Here we employ these manifestations of coherence to determine when mesoscopic quantum transport can be captured by a classical model based on stochastic jumps, and when such a model breaks down, implying nonclassical behavior. To this end, we focus on a minimal model of a double quantum dot coupled to two thermal reservoirs. In this system, quantum tunneling induces Rabi oscillations and results in both entanglement and nonlocality, as well as TUR and KUR violations. These effects, which describe the breakdown of a classical description, are accompanied by a peak in coherence. Our results provide guiding principles for the design of out-of-equilibrium devices that exhibit nonclassical behavior.


I. INTRODUCTION
The inability of classical physics to account for experimental observations on the atomic scale resulted in the emergence of quantum mechanics in the last century [1]. Quantum physics exhibits behavior that may be strikingly different from that of classical physics. One fundamentally distinguishing feature is quantum coherence, i.e., the ability to maintain a superposition of different quantum states. Coherence can be understood as a resource in quantum information processing [2], to realize tasks such as quantum computation or quantum metrology. Coherence also plays an important role in the thermodynamics of quantum systems and influences fluctuations of work and entropy production [3][4][5].
Despite the enormous interest in quantum advantages for technological applications, the difference between quantum and classical behavior has not been fully understood yet. One particular instance is that of quantum transport and thermodynamics in open quantum systems. Such systems may usually be described by multiple coupled sites, in contact with thermal reservoirs at different temperatures and/or chemical potentials (Fig. 1). After a long time elapses, this system will reach a nonequilibrium steady state (NESS), characterized by currents flowing between the reservoirs. This setup forms the basis for a variety of devices and applications, including thermoelectrics [6,7], autonomous engines [8], and nanoscale sensors [9]. * kacper.prech@unibas.ch which may have an important technological impact as it can be exploited for, e.g., quantum cryptography [25,26] or to provide advantages in sensing applications [27][28][29]. Entanglement, though, depends solely on the quantum state of the system at a given time; i.e., it is a static quantity.
The second manifestation is provided by violations of the so-called thermodynamic uncertainty relations (TURs) [30][31][32][33] and kinetic uncertainty relations (KURs) [33], as well as their recently proposed unification, the TKUR [34]. These are bounds on the signal-tonoise ratio of currents in classical systems. In quantumcoherent systems, these bounds no longer need to be respected, so that TUR or KUR violations can be interpreted as manifestations of quantum coherence. In contrast to entanglement, the TUR and KUR bounds depend on two-time correlations. These two types of manifestations of coherence are therefore complementary: while entanglement can be regarded as a static manifestation of coherence, TUR and KUR violations are dynamical in nature.
Here we focus on a minimal model, containing just two sites and two reservoirs (see Fig. 1). In this setup, both NESS entanglement [10][11][12] as well as TUR violations [16,17] are present. Investigating these manifestations of coherence together, and comparing to a stochastic jump model, allows us determine when quantum coherence results in a behavior that cannot be captured by a classical model. We find that the buildup of coherence due to quantum tunneling can become large enough to generate not only entanglement but even nonlocality, in contrast to previous studies [11]. Furthermore, the Rabi oscillations induced by quantum tunneling result in a reduction of fluctuations that may overcome both the TUR as well as the KUR. To the best of our knowledge, such KUR (and TKUR) violations have not been reported in quantum systems before. Interestingly, both manifestations of coherence occur when the tunnel coupling between the sites is comparable to the systembath coupling, where the coherence in the NESS shows a local maximum. In addition, entanglement is found for strong tunnel couplings, where coherence exhibits its global maximum. In this regime, however, the dynamics is well captured by a classical jump model that involves the entangled eigenstates of the two-site system. Thermokinetic uncertainty relations may therefore not be violated in this strong-coupling regime.
For concreteness, we consider a double quantum dot, where electrons can hop between two levels. The basic model parameters are outlined in Fig. 2. In order to contrast classical models and quantum-coherent models, we employ different quantum master equations [35,36] (local vs global) and investigate the effect of interactions. For a noninteracting system, we benchmark local and global quantum master equations with nonequilibrium Green's functions (NEGFs). We note that we benchmark not only average quantities, depending on the steady state, but also current fluctuations which depend FIG. 2. Sketch of the system. Two spin-polarized quantum dots (QDs) (ℓ = L, R) with on-site energies ϵ are tunnelcoupled with strength g and exhibit an interdot Coulomb interaction U . Each QD is weakly coupled to a fermionic reservoir with temperature T ℓ and chemical potential µ ℓ . The coupling between system and reservoir is denoted by γ ℓ . on two-time correlation functions. This further clarifies the regimes of validity of the local and global master equation which have been debated recently [35,[37][38][39].
In the presence of interactions, we employ a thermodynamically consistent semilocal master equation [36,40], as well as NEGFs. While our results are valid for the two-site model sketched in Figs. 1 and 2, the insights gained in the breakdown of a classical description carry over to more complicated transport scenarios, serving as guiding principles in designing out-of-equilibrium devices that can behave nonclassically. We illustrate this by investigating a chain of three quantum dots.
This paper is organized as follows. The system we consider and the models we employ for its description are discussed in Sec. II. The manifestations of coherence are introduced in Sec. III. In Sec. IV, we present the results for noninteracting electrons and in Sec. V we discuss the role of interactions. Section VI contains a discussion on more complicated transport scenarios, including a chain of three quantum dots. Conclusions are provided in Sec. VII.

A. Double quantum dot
We consider two quantum dots (QDs) connected in series, labeled ℓ = L, R [41][42][43]. Each QD is weakly coupled to a fermionic reservoir with temperature T ℓ and chemical potential µ ℓ , respectively (Fig. 2). As shown previously, this system exhibits both manifestations of coherence we are interested in: First, by driving an electron current through the system it is possible to generate an entangled state [10]. Second, fluctuations of the same current may be suppressed below the classical TUR bound [16]. We consider spin-polarized electrons, such that at most one electron may reside on each QD. The Hamiltonian of the system readŝ where ϵ is the onsite energy of the QDs, g the interdot tunnel coupling, U the Coulomb repulsion between electrons, andĉ ℓ (ĉ † ℓ ) the annihilation (creation) operator of the electron on QD ℓ, respectively. The effect of unequal onsite energies is considered in Appendix A.
The effect of the reservoirs is described by multiple methods, including local and global Lindblad master equations as well as NEGFs. We now introduce these methods for noninteracting electrons (U = 0); the interacting case is treated in Sec. V. The strength of the coupling between QD ℓ and the corresponding reservoir (assumed to be energy-independent) will be denoted by γ ℓ for all methods. Throughout, we consider a voltage bias that is applied symmetrically such that µ L = ϵ + eV /2 and µ R = ϵ − eV /2.

B. Lindblad master equations
The main method we use to investigate the dynamics of the double QD is provided by Lindblad master equations. These rely on Born-Markov approximations that for our system are justified when (we set k B = ℏ = 1 throughout) for any ℓ and ℓ ′ and for both signs. This condition ensures that the relevant bath properties are flat around the Bohr frequencies ϵ ± g. To ensure thermodynamic consistency, we follow a recently derived framework [36,40]. Both the local as well as the global master equation that we employ can be written in the form where L ℓ is the Lindblad superoperator that accounts for the dissipation to the bath ℓ. For the local master equation, the Lindblad superoperators are given by [36] where and n ℓ := n ℓ (ϵ), with the Fermi-Dirac distribution When the Born-Markov approximations are justified, the local master equation (ME) may be employed when [36] g ≪ max{T ℓ , |ϵ − µ ℓ |} (local ME).
In analogy to Eq. (2), this condition ensures that the relevant bath properties are flat across both Bohr frequencies. In this work, we use the local master equation extensively, as its regime of validity coincides with the regime where both manifestations of coherence appear.
In addition, we make use of the global master equation which relies on the secular approximation [44]. The Lindblad superoperators in the global approach read [36] where and n s ℓ := n ℓ (ϵ s ). The annihilation operatorsĉ ± correspond to the eigenmodes of Hamiltonian (1). The secular approximation is justified when [44] g ≫ max{γ L , γ R } (global ME).
In both Lindblad master equations, the steady-state average current flowing across the system can be computed with the expression is the operator of the number of electrons in the double quantum dot. To determine the current fluctuations, we resort to the method of full counting statistics discussed in Appendix B. We note that in both master equations, we neglect a Lamb shift that results in a shift of the on-site energies.

C. Nonequilibrium Green's functions
For noninteracting electrons, all relevant quantities can be solved exactly using the method of NEGFs. The transmission function T (ω), which specifies the transition rate of electrons with energy ω, is a fundamental quantity to investigate transport across the system. For our system it is given by [17,45] where we assumed the wideband approximation [46], i.e., energy independence of the coupling strength γ ℓ and omitting any Lamb shift. The transmission function allows us to compute the average of the current and its fluctuations, given by [47][48][49] and The elements of the density matrix of the system can be computed from the lesser Green's functions together with Wick's theorem. More information on the NEGF formalism can be found in Appendix A 3.

D. Classical model
The last method that we employ is a classical, stochastic model for the double QD system [15,50]. We consider the time evolution of the vector of probabilities ⃗ p = [p 0 , p L , p R , p D ] T , where entries denote, respectively, the probability that the system is empty, only the left QD is occupied, only the right QD is occupied, and both QDs are occupied. In the classical model, the time evolution of ⃗ p is governed by the rate equation where the elements W i̸ =j of the 4 × 4 matrix W are the transition rates from state j to state i, and W ii = − j̸ =i W ji . The transition rates between states with a different total number of electrons are taken from the local master equation, i.e., W L0 = W DR = n L γ L , W 0L = W RD = (1−n L )γ L , and similarly for L ↔ R. The transition rate describing hopping between the dots can be obtained from perturbation theory (cf. Appendix A 4) and reads The average current can be evaluated as To compute the fluctuations of the current we resort to full counting statistics (cf. Appendix B), as in the case of the Lindblad master equations.

A. Coherence
We study coherence in the quantum stateρ of the system with respect to the basis of the particle occupation number of each quantum dot. A local Fock basis {|n L , n R ⟩ := (ĉ † L ) n L (ĉ † R ) n R |0, 0⟩} is constructed such that the first (second) number denotes the occupation of the left (right) QD. Then, the density matrix of the system can be written in the basis {|0, 0⟩, |1, 0⟩, |0, 1⟩, |1, 1⟩} as a matrixρ It has a single off-diagonal element α, corresponding to a single electron being in a superposition between the left and the right dots. Any other off-diagonal element is strictly zero due to the charge superselection rule [51,52]. Throughout, we will quantify coherence with the absolute value |α|, which is equivalent to one half of the l 1 -norm of coherence [53].

B. Entanglement
The system sketched in Fig. 2 and its modifications have attracted a significant amount of attention because they may feature entanglement in the steady state [10][11][12][13][14]. This entanglement may arise from two different mechanisms. For tunnel couplings comparable to the system-bath coupling, the entanglement is generated by a current passing through the double dot [12]. This scheme to produce nonseparable states can be extended to multiqubit machines [54]. It can also be improved by local filtering [55], and it can be mediated by a cavity mode [56].
For tunnel couplings much larger than the temperature, only the lowest energy eigenstate is occupied. This state is an entangled state such that the thermal state becomes entangled. In this paper, we are mainly interested in the generation of entanglement by passing a current through the system because this is the regime where we also observe other manifestations of coherence.
We emphasize that the introduced scheme generates mode entanglement between two QDs, where the number of particles plays the role of the degree of freedom [52,[57][58][59]. This is qualitatively different from entanglement between two particles, such as spin entanglement in electronic systems. We note that the entanglement of a single electron [60] has been debated in the literature due to the charge superselection rule. However, following [58], having two copies of the system allows to overcome this issue and makes the entanglement of a single electron useful.
To quantify the amount of entanglement we use the entanglement measure of concurrence [61,62]. For the present system, the concurrence is given by From this expression, we can already anticipate that Coulomb interactions have a favorable effect on the entanglement, because they reduce p D . This has a twofold effect as it allows for a larger occupation in the singleelectron subspace, which may result in a larger α.

C. Thermodynamic and kinematic uncertainty relations
The TUR and KUR are two different inequalities that bound the signal-to-noise ratio of a current in classical systems. The TUR has been developed in the field of classical stochastic thermodynamics [63][64][65][66][67], which is a framework to describe fluctuations of thermodynamic quantities, such as work or heat, in nanoscale systems. It provides the bound in terms of the average entropy production rate, which quantifies the irreversibility of a process. The TUR can be applied to obtain a trade-off between power production, constancy, and efficiency in stochastic heat engines [68] and to infer the efficiency of molecular motors [69].
The KUR was derived using the fluctuation-response inequality for out-of-equilibrium stochastic dynamics [70] and finds applications outside of the thermodynamic settings. In the case of the KUR, the signal-to-noise ratio is bounded by the average dynamical activity, which is a measure of the total rate of transitions between the states of the system. The question whether the KUR can be violated in a quantum system has, to the best of our knowledge, not been explored. We note that in contrast to entropy production, which can be derived using any model, subtleties arise when computing the dynamical activity that features in the KUR with different models. These are discussed below.
While the TUR provides a tight bound when dissipation is small, the KUR is more relevant for systems that are far from equilibrium. Either of these two different inequalities may thus provide a stronger bound on the fluctuations depending on the situation. As discussed below, these inequalities have also recently been combined into the TKUR [34].
The TUR quantifies a trade-off between an average of a current, ⟨I⟩, its fluctuations ⟨⟨I 2 ⟩⟩, and the total entropy production rate ⟨σ⟩. It was initially proven in the long-time limit for nonequilibrium systems that follow a time-homogeneous Markovian dynamics and obey local detailed balance. It is given by [30,31] In addition, several modified versions of the TUR have been derived, which include a finite time generalization [71,72], measurement and feedback scenarios [73], as well as time-dependent periodic [74] and nonperiodic [75] driving. In this work, however, we restrict our attention to the TUR given in Eq. (21). In the steady-state regime, entropy production comes solely from the dissipation of heat within the reservoirs [6], where J ℓ is the steady state heat current from the reservoir ℓ. How these heat currents can be computed from the different models is discussed in Appendix A. For reservoirs with equal temperature T , the entropy production simplifies to Moreover, Eq. 21 reduces to which only depends on the Fano factor ⟨⟨I 2 ⟩⟩/⟨I⟩. The KUR provides an alternative bound on the signalto-noise ratio [33], where ⟨A⟩ is the dynamical activity that quantifies the average rate of jumps the system undergoes.
It quantifies the average total rate of jumps between all states of the system. In the quantum model, the process of electrons moving between the left and right QDs is modeled by a coherent evolution, not classical jumps. For this reason, the dynamical activity for Markovian quantum master equations is usually defined such that it only quantifies the rate of jumps induced by Lindblad jump operators [76,77]. For the local master equation, this would read resulting in a strictly smaller dynamical activity. Employing ⟨A q ⟩ in the KUR could result in the following problem: In a regime where the classical model and the master equation provide the same probabilities and fluctuations (i.e., the system essentially behaves classically), we could get a violation of the KUR simply because ⟨A q ⟩ does not take into account the transitions between the left and right dots. To avoid such problems, we employ Eq. (26) throughout for the dynamical activity, where the rates W ij are always taken from the classical model (16), and the probabilities p j are obtained with the respective method that is employed. Finally, we consider the TKUR that combines the above bounds as [34] where f (x) is the inverse function of x tanh x. This inequality optimizes between the TUR and the KUR, and thus gives a stronger bound than both of them. When investigating violations of uncertainty relations, we use the quantifier where j = T, K, T K correspond to TUR, KUR, and TKUR, respectively. We refer to this quantity as j violation. If it is positive, the respective uncertainty relation is violated. We note that numerically comparing the values of these quantifiers between different types of violations is not meaningful.

IV. NONINTERACTING ELECTRONS
In this section, we provide our results for noninteracting electrons. In this case, there exist exact solutions obtained with NEGFs (cf. Appendix A 3), which serve as a benchmark. Before illustrating the manifestations of coherence, we consider the influence of the tunneling strength g on coherence, the current, and its fluctuations. This allows for illustrating the coherent nature of the dynamics and to compare the different models. All analytical expressions provided in this section are derived with the local Lindblad master equation. For analytical results obtained with the global approach we refer the reader to Appendix A 2.

A. Coherent dynamics: Comparing the models
We first discuss how the coherence in the system changes with the interdot tunnel coupling in the steady state. We find Interestingly, this expression exhibits a peak at g ≃ γ (with γ L = γ R = γ) [see Fig. 3(a)]. This peak in coherence can be understood as follows: The coherent tunneling induces Rabi oscillations between the left and the right dots, which may result in a buildup of coherence. For g ≪ γ, the decoherence induced by the bath suppresses any buildup of coherence akin to the Zeno effect. For g ≫ γ, the Rabi oscillations become very fast, suppressing coherence by phase averaging. For intermediate g, coherence can build up, resulting in a peak. We note that the location of this peak does not depend on temperature or chemical potential (which only enter in the Fermi-Dirac distributions n ℓ ). While the peak in coherence is captured well by the local master equation, this approach breaks down when g becomes of the order of T [cf. Eq. (7)]. In this regime, where the global master equation is justified, coherence grows again. As discussed above, this is a result of the entangled ground state being more populated than the excited state. In the large g limit, the steady state reduces to the ground state, a pure singlet state with maximal coherence. As illustrated in Fig. 3(a), we find excellent agreement between the exact NEGF solution and the Lindblad master equation solutions in their respective regime of validity.
Additionally, the current and its fluctuations provide insight into the effect of coherence on the dynamics. For the current, the local master equation yields The average current is illustrated in Fig. 3(b). Starting from small g, its value grows as g 2 and saturates when the tunnel coupling g becomes larger than the systembath coupling γ, which becomes the bottleneck for transport. Similarly to the coherence, when g is of the order of T we observe a breakdown of the local approach. For larger values of g, where the global approach describes the dynamics well, the average current decreases, as the eigenenergies of the Hamiltonian leave the bias window, and transport of electrons is impeded. In addition, the plot includes the average current predicted by the classical model (cf. Appendix A 4). It is in complete agreement with the local master equation, and as such is accurate in the regime of validity of the local approach.
The current fluctuations are given by which is illustrated in Fig. 3(c). Comparison with the NEGF solution shows that the local and global master equations reproduce the exact current fluctuations in their respective regime of validity. To the best of our knowledge, this is the first time that these regimes of validity have been benchmarked by quantities that do not only depend on the quantum state alone but instead depend on two-time correlation functions. As a function of g/γ, the fluctuations display analogous features to the average. However, when g is of the same order as γ, the fluctuations are suppressed below the value of the classical model (see also Ref. [15] for a similar observation). The regime characterized by the peak in coherence thus coincides with the regime where the fluctuations of the current are suppressed relative to the classical model. These features underlie the manifestations of coherence that we focus on below.
All three manifestations are presented in Fig. 4 as a function of g and eV . The range of values of g coincides with the window of the tunnel couplings that exhibit both a peak of coherence and reduction of current fluctuations. This is directly demonstrated in Figs. 4(b) and 4(c), where we compare |α| with C, V T , V K , and V T K on two cross sections of Fig. 4(a). Both plots show excellent agreement between the local master equation (solid color lines) and NEGFs (corresponding dashed black lines). The regime where g and γ are comparable thus features both entanglement, a static manifestation of coherence describing the steady state, as well as TUR and KUR violations, dynamical manifestations of coherence which arise from the nonclassical dynamics of the system. Crucially, we see how these two manifestations of coherence do not always appear together. As discussed above, in the regime of strong g, entanglement is found. However, the dynamics of the system are well described by the global master equation, a classical rate equation involving the entangled eigenstates of the Hamiltonian. Consequently, we do not find TUR and KUR violations.
There are also regimes where entanglement and TUR violations are found to partially overlap, such that a crossover can be achieved by increasing the voltage bias. Specifically, the TUR is only violated for small voltage bias, whereas entanglement requires large voltage bias. Increasing eV gives rise to increasing entropy production, which causes the TUR to become less tight. Entanglement grows with increasing current, which itself increases with the voltage. For large voltages, we also observe KUR violations due to the suppression of current fluctuations. In agreement with previous findings V K ≈ 0.216, TUR violations are observed close to equilibrium while KUR violations appear far from equilibrium. As the TKUR is tighter than both the TUR and the KUR, its violations encompass both the TUR as well as the KUR violations.
We end this section by considering the usefulness of the entangled quantum state to achieve nonlocality [24] and to perform quantum teleportation [78] (see Appendix D). This has been done previously in a similar system [11]. For noninteracting electrons, we find a teleportation fidelity of f = (7 + √ 5)/12 ≈ 0.77 (cf. Appendix D 2), which is above the classical limit of 2/3. As expected, this corresponds to the value found in Ref. [11] in the noninteracting case. Concerning nonlocality, in our system with noninteracting electrons we do not find any violation of the Clauser-Horne-Shimony-Holt (CHSH) inequality [79] (cf. Appendix D 1), in agreement with the findings of Ref. [11]. As discussed in Sec. V C, inclusion of Coulomb interactions allows for a higher teleportation fidelity and even CHSH violations, implying Bell nonlocality. This is in contrast to the results of Ref. [11], where no chemical potential was present and population inversion was obtained using negative temperatures. Coulomb interactions between the electrons change the dynamics of the double QD system by incorporating an additional energetic cost for occupying both QDs. In this section we show how the presence of interactions affects the manifestations of coherence.

A. Methods
We investigate the dynamics of the system using methods analogous to the noninteracting case, adjusted to take into account the effect of Coulomb interactions between electrons.

Master equation
In the previous section we have identified the regime of interest to be where the peak of coherence appears together with all the manifestations of coherence. Since the same holds for interacting electrons, we restrict our attention to the regime where tunnel couplings g are comparable to γ ℓ . A Lindblad master equation suitable for treating this regime in the presence of Coulomb interactions (termed semilocal) has recently been microscopically derived [36]. The corresponding superoperators read where andl ̸ = ℓ. It is valid in the regime The key difference between the semilocal master equation and the local one, which we used for the noninteracting case, is that now the jump rates from each QD to the adjacent reservoir depend on the occupation of the other QD. Importantly, the results from the semilocal master equation reduce to the results from the local master equation in the limit U → 0 [36].
In order to compute the average current and its fluctuations we proceed analogously to the case of noninteracting particles. The average current in the steady state is given by Eq. (11), and the fluctuations are obtained with full counting statistics (Appendix B).

Classical model
As for noninteracting electrons, we derived a classical model of the form of Eq. (16) using the perturbative approach described in Appendix A 4. The matrix elements describing jumps of electrons between the reservoirs and the system read W L0 = γ L n L , W 0L = γ L (1 − n L ), W DR = γ L n U L , W RD = γ L (1 − n U L ), and similarly for L ↔ R. These are the same rates that appear in the semilocal master equation [cf. Eq. (35)]. The interdot tunnel rate is given by and W RL = W LR .

NEGFs
The effect of Coulomb interactions in the NEGF can be introduced via a many-body self-energy. This is an inprinciple exact procedure. In practice, one must, however resort to approximate self-energies, calculated, e.g., via many-body approximations. Here, we consider a rather simple approximation, namely the second Born perturbative scheme [80][81][82] (for details we refer the reader to Appendix C 2). The range of U/g where the second Born approximation gives a reliable description of electronic correlations is limited and much smaller than the range we consider with the Lindblad master equation in the remainder of the paper. Nonetheless, it may still be useful to compare the two methods in a restricted domain where both approaches produce meaningful results.
We thus notice that in this limited U/g range the two treatments give results that are in good mutual agreement for the coherences, while discernible discrepancies occur for the concurrence, due to well-known possible shortcomings of the second Born approximation when determining double occupancies (see Appendix C 2).

B. Enhanced manifestations of coherence
Since Coulomb interactions introduce an energetic cost of having both QDs simultaneously filled, the probability of double occupation, p D , is suppressed and vanishes in the limit U → ∞. This reduced occupation of the doubly filled subspace generally has a beneficial effect on the amount of coherence, as it increases the probability for the subspace with a single electron on the dot. However, as explained below, Coulomb interactions do not always result in an increased amount of coherence. The expression for coherence reads where we have introduced and ∆ l = n l − n U l . In the limit U → ∞ the coherence reduces to .
(42) As illustrated in Fig. 5(a), increasing U results in a larger peak in coherence that is shifted to smaller values of g/γ. From the definition of concurrence, given by Eq. (20), we anticipate that Coulomb interactions increase the amount of entanglement in the system for two reasons: First, because p D is smaller, and second, because this may increase the off-diagonal element α. This is illustrated in Figs. 5(b) and 5(c), where we survey the influence of U/g on the coherence and concurrence, respectively. In Fig. 5(b) coherence shows moderate growth with the strength of the Coulomb interaction. For comparison, we feature the results of the semilocal master equation (black) and of NEGFs (red) on the restricted interval, demonstrating a good agreement between the two methods. A similar plot showing the concurrence is presented in Fig. 5(c). For concurrence we observe a strong impact of Coulomb interactions, which indicates that the growth of concurrence comes predominantly from the suppression of p D . We note that the agreement between the semilocal master equation and NEGFs is better for coherence than for concurrence. We attribute this to the fact that the concurrence relies not only on single particle quantities but requires the computation of p D . As discussed in Appendix C 2, this quantity may be more sensitive to the employed approximations.
In addition to entanglement, Coulomb interactions impact the dynamical manifestations, TUR and KUR violations, which is illustrated in Fig. 6. In the series of panels we show how C, V T , and V K evolve as U/T increases. Concurrence and TUR violations are considerably enhanced by interactions, while the enhancement of KUR violations is less prominent. While the range of g/γ where entanglement is present is extended, the same is not true for the remaining manifestations. This is a consequence of the fact that larger C is caused predominantly by the reduction of p D as opposed to an increase in coherence. In addition, the values of g/γ that maximize V T and V K become smaller and align with the optimal value for C. This is similar to the observed shift in the peak of coherence. The main drawback of looking at the entanglement solely through the prism of any entanglement measure, such as concurrence, is that it does not provide a definite measure of the usefulness to achieve nonlocality [24] or perform quantum information tasks, such as quantum teleportation [78]. The problem of generating operational nonclassicality in the class of systems akin to the double quantum dot was systematically analyzed by Bohr Brask et al. [11], where no nonlocality was found. In their approach, however, the system was brought out of equilibrium via a temperature gradient alone but allowing for negative temperatures to describe population inversion. Here we show that, by considering a bias voltage, Coulomb interactions may result in nonlocality.
The optimal regime to maximize concurrence and observe nonlocality is provided by the following conditions: First, jumps from the right reservoir into the system are completely suppressed, i.e., n R = n U R = 0. Second, jumps from the left reservoir are suppressed only when there is already an electron in the system, i.e., n U L = 0 and n L = 1. This is achieved in the limits U → ∞ and eV → ∞, with U/eV → ∞.
Under these conditions, the probability of the double occupation p D vanishes and the coherence reduces to The expression for concurrence simplifies to Upon examining ∂C ∂g , we find the maximal value C = √ 2/2 ≈ 0.71 when γ R /g = 2 √ 2 and γ L /g → ∞. This is significantly larger than ( √ 5 − 1)/4 ≈ 0.31, which we found to be the maximum in the noninteracting case. If the couplings to the baths are assumed to be equal, the largest concurrence occurs at γ/g = 2 √ 3, resulting in C = √ 3/3 ≈ 0.58. We note that in the absence of a chemical potential (even allowing for negative temperatures), Coulomb interactions do not alter the maximal value for the concurrence [11].
Motivated by the substantial increase in concurrence due to Coulomb interactions, we revisit the question of Bell nonlocality (cf. Appendix D 1), focusing on the CHSH scenario [79], which was also considered in Ref. [11]. In this scenario, if the statistics of the outcomes of measurements on a quantum state admit a local hidden-variable model, then the CHSH quantity obeys CHSH ≤ 2 [24]. Any violation of this bound is a signature of nonlocality, and the maximal CHSH = 2 √ 2 can be achieved with maximally entangled Bell states. For the double quantum dot, in the regime discussed above, we find the maximal value CHSH = 2 3/2, which is achieved for the same set of parameters γ L , γ R , and g that saturate concurrence. When γ R /g = 2 √ 2, CHSH displays nonlocality for γ L /g > 5. 66.

VI. MORE GENERAL MODELS
Up to this point we have been concerned with the minimal model consisting of only two sites, which allows to clearly identify the regimes of nonclassical behavior. The fact that similar conclusions hold for both interacting and noninteracting electrons implies that our results may serve to anticipate the behavior of more general mesoscopic systems. To illustrate this, we investigate TUR violations in a chain of three QDs, where electrons in the system may coherently traverse between the neighboring QDs. The corresponding Hamiltonian is given bŷ where the subscript M refers to the middle QD, and g 1 (g 2 ) denotes the coherent coupling between the left and the middle (the right and the middle) QDs. We investigate the chain with the NEGFs using analogous methods as those presented in Appendix A 3. As anticipated from our simple model, we find nonclassical behavior where both coherent couplings, g 1 and g 2 , are of the order of the coupling to the baths γ. This is illustrated in Figs. 7(a) and 7(b), which show TUR violation V T and the amount of coherence in the system, respectively. As a quantifier of the amount of coherence we use one-half of the ℓ 1 -norm of coherence, here denoted by C ℓ1 , which reduces to |α| in the two-site model. The peaks of V T or C ℓ1 are visualized in Fig. 7(c), where we show the diagonal cross-sections with equal coherent couplings, g 1 = g 2 = g, of Figs. 7(a) (black) and 7(b) (red), respectively. These graphs further demonstrate that the TUR violation and coherence are maximized for g 1 = g 2 both being comparable to γ. Moreover, in the chain of three QDs we observe stronger TUR violation than in a double quantum dot (DQD) operating under the same set of parameters (see Fig. 4). This indicates that nonclassical behavior is more easily accessible in longer chains. The fact that our conclusions, drawn from the investigation of the DQD, hold in the chain of three QDs provides evidence for the applicability of these results in more general systems. We note that these results are also consistent with Ref. [84,85], where TUR violations in longer chains of QDs were investigated.

VII. CONCLUSIONS AND OUTLOOK
In this paper we investigated the nonclassical behavior of a mesoscopic open quantum system by analyzing two complementary manifestations of coherence: Entanglement, a static manifestation, and violations of thermokinetic uncertainty relations, dynamical manifestations of coherence. We considered a double quantum dot weakly coupled to fermionic reservoirs, where the transport of electrons between the dots can be modeled either by coherent tunneling, or by stochastic jumps.
We identified the regime where the system exhibits truly nonclassical behavior: On the one hand, coherence exhibits a peak, which may result in an entangled state and, in the presence of interactions, even nonlocality. On the other hand, coherent tunneling suppresses the fluctuations of the current below the classical model, allowing for violations of the TUR and the KUR. This behavior occurs in the regime where the coherent tunnel-coupling is of the order of the coupling to the baths, which is captured by the local master equation. This is different from the regime of large tunnel couplings, which is captured by the global master equation, where the state of the system may be strongly entangled, but the dynamics of the system can be expressed through a classical rate equation, and no TUR or KUR violations are present.
Our systematic characterization of the electron transport across the system performed with different master equations and NEGFs illustrates the need to go beyond the steady state in order to fully capture the nonclassical behavior in mesoscopic transport. While the steady state captures the entanglement, the violations of thermokinetic uncertainty relations cannot be captured by the steady state alone. Those different manifestations of coherence are not equivalent and do not need to coincide. We note that in the regime of nonclassical behavior, not only the quantum state of the DQD but also the current fluctuations, which depend on two-time correlations, are well captured by the local master equation. While our results are obtained for a specific model, we believe that these conclusions are relevant for a broad range of open quantum systems and provide guiding principles for the design of out-of-equilibrium devices that exhibit nonclassical transport. This is evidenced by the fact that our predictions about the nonclassical behavior, drawn from the noninteracting two-site model, also hold for interacting electrons and in the chain of three QDs.
Our results may be probed experimentally in a number of different platforms, including electronic QDs, which can be implemented, e.g., using two-dimensional electron gases [86,87] or nanowires [88,89], as well as superconducting qubits [90][91][92]. While measurements of the electric current and its fluctuations across a DQD is a standard experimental technique [93], accessing the quantum state of the DQD is more challenging. Alternatively, the density matrix and, therefore, the entanglement may be probed in an implementation of our model using two coupled superconducting qubits [90]. In the main text of the paper we restrict our attention to the case with the same onsite energies of both QDs. In this appendix, however, we provide a derivation of the results in the general case of noninteracting electrons. The Hamiltonian operator of the system is given bŷ where ϵ ℓ is the on-site energy of the QD ℓ. In the noninteracting case, the coefficients of density matrix (19) can be expressed as where we have used a corollary of Wick's theorem [94],

Local master equation
In the local master equation satisfying thermodynamic consistency the Lindblad superoperators appearing in Eq. (3) are given by [36] where now n ℓ := n ℓ (ε), andε = (ϵ L + ϵ R )/2. In order to find the density matrix in the steady state, we consider the vector where Using Eq. (A2) we obtain Using the formula for the concurrence, given by Eq. (20), we obtain the following expression, if this quantity is positive and zero otherwise.

a. Average current
In the steady state, the average current passing through the system from the left reservoir can be obtained from the expression given in Eq. (11). Equivalently, one may compute the current as an expectation value of the current operator, In the steady state, we find

. Heat current
We define the heat currents in the thermodynamically consistent framework [36,40], which guarantees the ther-modynamically consistent heat flow and maintains the accuracy of the local form of the master equation. In this framework the average heat current from the reservoir ℓ to the double quantum dot is given by [36] where we have introduced a thermodynamic bookkeeping HamiltonianĤ We find and We note that if ϵ L = ϵ R , which is assumed in the main text of this paper, the steady-state heat currents reduce to the ones obtained from the conventional definition, where in Eq. (A13) the HamiltonianĤ is used instead of H T D .

c. Current fluctuations
In order to investigate the fluctuations of the current, we resort to the method of full counting statistics. All the details about this procedure are explained in Appendix B.
The counting-field-dependent Liouvillian L(χ), introduced in Eq. (B5), is obtained by inserting the terms e ±iχ into the Lindblad superoperator [Eq. (A4)] of the local master equation, where L is the Liouvillian in the absence of a counting field given in Eq. (3). Writing L(χ) in the basis The expression for the current fluctuations is given by
From the global master equation, we find the equations of motion The corresponding steady-state solution is given by and the remaining expectation values vanish, i.e., ⟨ĉ † +ĉ− ⟩ = ⟨ĉ † −ĉ+ ⟩ = 0. The density matrix elements in the local basis can be obtained from Eq. (A2), using the relations The coherence then reads whereas the concurrence is given by

a. Average current
As for the local master equation, the average may be obtained from Eq. (11). We note that for the global master equation, we may not use a current operator [35]. We find The heat current from reservoir ℓ may be written as [36] resulting in wherel ̸ = ℓ.

c. Current fluctuations
Similarly to the local master equation, we find a χdependent Liouvillian: In the global master equation, the natural choice of the basis is given by {p 0 , p + , p − , p D }, where p ± = ⟨0, 0|ĉ ±ρĉ † ± |0, 0⟩. In this basis the Liouvillian reads (A35) Using the method of full counting statistics outlined in Appendix B, we find

Nonequilibrium Green's functions
The key building blocks in the steady-state description with NEGFs are the retarded G r (ω) and the lesser G < (ω) Green's functions [95,96]: We use boldface symbols to represent 2 × 2 matrices in the site indices (L, R) of the QD system. 1 is the identity matrix, the Hamiltonian of the system [cf. Eq. (A1)] is represented by and G a = (G r ) † . In the absence of Coulomb interaction the retarded or lesser self-energy Σ r/< (ω) describes the coupling between the double quantum dot and the reservoirs. In the wideband limit we have [46,96,97] and From Eqs. (A37) and (A38) we find the following four elements of G < (ω), given by The expectation values ⟨ĉ † ℓĉℓ ′ ⟩ (ℓ, ℓ ′ = L, R) are obtained by taking the inverse Fourier transform of the corresponding lesser Green's functions [95,96] where we use the convention Cumulants of the electron current through the system can be obtained from the transmission function T (ω). The expression for the average current reads [47][48][49] and the fluctuations are given by The transmission function itself can be expressed with the Green's functions [96] T where For the double quantum dot system it is given by [17,45] In the same framework, the heat current from reservoir ℓ reads wherel ̸ = ℓ.

Stochastic model of the system
To motivate the classical model, given by the rate equation in Eq. (16), we consider the local master equation with L given in Eqs. (3) and (4). We now introduce Nakajima-Zwanzig superoperators P and Q = I − P, such that P projects onto the population subspace ofρ. Here I is the identity superoperator. We note that P and Q are projectors, i.e., P 2 = P, Q 2 = Q, and QP = 0. We also decompose the Liouvillian of the system into L = L 0 + V, where V is the contribution due to the inter-dot tunneling, i.e., Vρ = −ig[ĉ † Lĉ R +ĉ † Rĉ L ,ρ]. For small g we can treat V as a perturbation. In the following derivation, we will use the relations By acting with P and Q on Eq. (A51) we obtain the equations and The formal solution to the second equation reads By inserting L = L 0 + V, applying Eq. (A52), and using the properties of Nakajima-Zwanzig projectors, we obtain the formal solution where we have also made a substitution s → t − s. The eigenvalues of V and L 0 are proportional to g and γ ℓ , respectively. For g ≪ γ ℓ , we can treat V as a perturbation and neglect its contribution to the exponential e (QL0Q+QVQ)s . This exponential thus decays on the timescale 1/γ ℓ . We now assume that the diagonal ele-ments of the density matrix remain approximately constant during this timescale, an assumption which is not justified (see below for more information). This allows us to substitute Pρ(t) for Pρ(t − s) in Eq. (A58) and to extend the integral to infinity. Equation (A58) then reduces to a Markovian master equation for Pρ, where is the Drazin inverse of L 0 .
Equation (A59) can be cast into the form given in Eq. (16) with the stochastic matrix All entries in W are identical to the ones in the local master equation [cf. (A18)] except the element W LR = W RL , which comes from the term PVQL −1 0 QVP (here we assumed δ = 0 for simplicity).
We stress that the approximation giving rise to Eq. (A59) is not justified. It assumes that the diagonal elements remain constant over the timescale 1/γ ℓ . This is, however, exactly the timescale over which the diagonal elements do change. For our purposes, this is not problematic since our aim is to derive an analogous classical description for the system. We do not claim that this classical analog is exact in any limit.
Interestingly, we find the steady-state populations in the classical model to be the same as in the local master equation given in Eq. (A7) for all choices of parameters. The same holds for the average current, which agrees with expression (31) obtained with the local master equation. The fact that the derived rate equation reproduces the average behavior of the local master equation further justifies its use as an analogous classical model.
To obtain the current fluctuations, counting fields may be introduced in analogy to Eq. (A18). Using full counting statistics, we then obtain (A62)

Appendix B: Full counting statistics
In order to obtain the expressions for the current fluctuations, we employ the technique of full counting statistics. In this framework the density matrix of the system is resolved intoρ(n, t), the (unnormalized) density matrices given that a net number of n electrons tunneled from the left reservoir into the DQD during the time interval [0, t] [98][99][100]. The probability that n electrons have tunneled can then be obtained by The original density matrix can be recovered by taking the sumρ The cumulants of the net number of transported electrons can be obtained from the cumulant generating function S(χ, t), given by and χ is known as the counting field.
The time evolution ofρ(χ, t) can be found by inserting the χ-dependent term, which tracks electron jumps, into the original Lindblad master equation [100]. This results in a counting-field-dependent master equation, The specific form of L(χ) depends on the choice of the master equation, and the original Liouvillian can be recovered with L = L(0). The formal solution to Eq. (B5) is given byρ The cumulants of the electron current in the steady state can be obtained from the scaled cumulant generating function of n in the long-time limit, where λ(χ) is the eigenvalue of L(χ) with the largest real part. The average of the current, ⟨I⟩ = ⟨⟨I⟩⟩, and its fluctuations, ⟨⟨I 2 ⟩⟩, can be evaluated by taking k = 1, 2, respectively.
Since a direct of computation of λ(χ) is not straightforward for our system, we make use of the scheme described in [101]. We express the Liouvillian L(χ) as a matrix L(χ) in a basis that can be freely chosen and may be different depending on whether we use the local or global master equation. Let P χ (x) be the characteristic polynomial where I is the identity matrix. Upon expanding in x we obtain where Here M is the dimension of the matrix L(χ). By definition, the eigenvalue λ(χ) of L(χ) (with the largest real part) fulfills It follows that In order to obtain the right-hand side, we first expanded P χ (λ(χ)) using Eq. (B9), used λ(0) = 0, as well as Eq. (B7). We thus find The second cumulant of the current may be obtained in a similar fashion from resulting in (B15) Appendix C: Interacting electrons When investigating the system interacting electrons we restrict ourselves to the case with equal onsite energies of quantum dots, where the Hamiltonian of the system is given by Eq. (1). In the presence of Coulomb repulsion, the coefficients of the density matrix in the basis already introduced in Eq. (A2) are given by (C1)
The steady state of the system is given by the eigenvector corresponding to the zero eigenvalue of the Liouvillian. TheLiouvillian can be computed from Eq. (C12) evaluated at zero counting field χ = 0. We find and ∆ l = n l − n U l . Expressions for coherence |α| as well as the concurrence C, which is defined in Eq. (20), can be obtained from the coefficients shown above.

a. Average electron current
For the average electron current through the system, given by Eq. (11), we find . (C7)

b. Heat currents
The heat current from bath ℓ is defined in complete analogy to the local master equation, i.e., Eq. (A13) [36]. In the presence of Coulomb interactions, the bookkeeping Hamiltonian readŝ resulting in and respectively.

c. Current fluctuations
The current fluctuations are obtained in complete analogy to the method used in the noninteracting case, i.e., full counting statistics (cf. Appendix B). The χ-dependent Liouvillian reads In the basis {p 0 , p L , p R , p D , Re[α], Im[α]} it is given by where γ U = − l 1−n l +n U l 2 γ l . We thereby find the current fluctuations (C13)

Nonequilibrium Green's functions
To determine the densities of the system via NEGFs, we can still proceed as in Appendix A 3. However, we have to account, within NEGF, for the Coulomb interaction via a many-body contribution to the self-energy [80][81][82]. As before, the key building blocks in the steadystate description with NEGF in the presence of interactions are the retarded G r (ω) and lesser G < (ω) Green's functions, defined in Eqs. (A37) and (A38) and for convenience repeated here: The self-energy Σ now has two parts, coming from the embedding (emb) and many-body (MB) contributions: where the embedding term has been discussed before in Appendix. A 3. The correlation effects due to interactions are introduced via the many-body self-energy Σ r/< For Σ MB we use the second Born approximation [97,102,103], keeping all Feynman diagrams up to second order. These equations are then solved self-consistently, using the Pulay scheme [104] to improve convergence. Selfconsistency ensures that general conservation laws [80][81][82] are obeyed (e.g., if I L is the current from lead L, then I R = −I L ). The embedding self-energy is independent of the self-consistency cycle and is calculated separately only once. For the particle density in the device region and the average current ⟨I⟩ coming from the left reservoir, the expressions are [96], respectively, Tr Γ L G < (ω) − 2πin L (ω)A(ω) .

(C17)
Here 2πA(ω) = i(G r (ω) − G a (ω)). Finally, the concurrence C is computed starting from the expressions for ⟨N ℓ ⟩ (ℓ = L, R), directly obtained from Eq. (C16), and from ⟨N LNR ⟩. The latter quantity is obtained according to a prescription introduced in Ref. [105]. We start with the defining expression for Σ. This in standard NEGF notation on the Keldysh contour reads in terms of the two-particle Green's function G 2 : Here and in Eq. (C18), 1 ≡ (ℓ 1 , z 1 ) (and similarly for 2, 3, 4), with ℓ 1 a site index and z 1 a time variable on the Keldysh contour.
We now make the site indices explicit in Eq. (C18), i.e., 1 → (ℓ, z), 2 → (ℓ ′ , z ′ ), 3 → (l,z), and set (ℓ ′ , z ′ ) = (ℓ, z + ) = (L, z + ). Furthermore, with nonlocal interactions, the interaction term becomes ν ℓl (z,z) = U (δ ℓL δl R + δ ℓR δl L )δ(z −z). Equation C18 then becomes l dz Σ Ll (z,z)Gl L (z, z + ) = −iU G LRLR (z, z + ; z + , z ++ ), where, here and in the following,l = L, R in the sums. Then, using the definition of G 2 in Eq. (C19), l dz Σ Ll (z,z)Gl L (z, z + ) and by applying the Langreth rules, we arrive at an expression for physical times: Finally, moving to frequency space, we have For computational convenience we separate the selfenergy into a correlation part and one coming from the Hartree-Fock interaction: (Σ < c,Ll (ω)G ā ℓL (ω) + Σ r c,Ll (ω)G < ℓL (ω)), where the subscript c on the self-energy indicates that it contains only the contribution from correlations. The applicability of this expression has been examined in Ref. [105]. There it was shown that some approximate but conserving many-body self-energies (among them, those obtained in the second Born approximation) do not always guarantee a non-negative value for local double occupancies ⟨N ℓ↑Nℓ↓ ⟩ in a spinful model with Hubbard interactions. Here, by contrast, we are considering a spinless model, with nonlocal interactions and correlations and in all applications of Eq. (C23) in the present paper, no negative correlations were encountered.
Appendix D: Operational nonclassicality Following Ref. [10], we consider the usefulness of the DQD steady state for Bell nonlocality and quantum teleportation.

Nonlocality
In the CHSH scenario [79] two players that share a bipartite quantum state, make measurements x and y, and obtain outcomes a and b. The conditional probability p(a, b|x, y) is nonlocal if it does not admit a local hiddenvariable model p(a, b|x, y) = λ p(λ)p(a|x, λ)p(b|y, λ) (D1) for any λ with an underlying distribution p(λ) [106]. In the CHSH scenario a, b, x, y ∈ {0, 1}, and p(a, b|x, y) can be explained by a local hidden-variable model whenever [79] CHSH = a,b,x,y (−1) a+b+xy p(a, b|x, y) ≤ 2.
For the density matrix of the double quantum dot, given in Eq. (19), the maximal value for the CHSH quantity is given by [11] For noninteracting electrons, we examine the maximum possible CHSH in the regime eV → ∞, where the corresponding Fermi distributions reduce to n L → 1 and n R → 0. Using the expressions for the coefficients of the density matrix, given in Eq. (A8), we do not find any parameters of the system that give rise to a violation of CHSH inequality.

Quantum teleportation
Quantum teleportation [78] is a protocol where an unknown quantum state is sent between two parties by the means of a shared entangled state and classical communication. In the original implementation, a maximally entangled Bell state |ψ − ⟩ = (|01⟩ − |10⟩)/ √ 2 was used, giving the average fidelity of teleportation f = 1. However, when the parties implementing the protocol do not share a Bell state, but any other less entangled stateρ, the fidelity of the teleportation is lower with the threshold value f = 2/3 for separable states. The average fidelity f = (1+2F )/3 can be computed from the so-called singlet fraction [83] where the optimization is performed over the unitary transformations U . For the density matrix of the double quantum dot, given in Eq. (19), we have that if its singlet fraction exceeds that of any separable state, i.e., F > 1/2, then it is given by [11] In the absence of Coulomb interactions F is saturated at the maximal bias, i.e., when eV → ∞. Using the expressions for the steady state of the density matrix, given by Eq. (A8), we find the corresponding singlet fraction, It is maximized for the symmetric coupling strength to the baths that fulfills g/γ = ( √ 5 − 1)/4, resulting in F = (3 + √ 5)/8 ≈ 0.65.
Under the conditions described in Sec. V C for interacting electrons, p D = 0, whereas |α| and p 0 are given by Eqs. (43) and (C5), respectively. The singlet fraction is then given by .