On optimal currents of indistinguishable particles

We establish a mathematically rigorous, general and quantitative framework to describe currents of non- (or weakly) interacting, indistinguishable particles driven far from equilibrium. We derive tight upper and lower bounds for the achievable fermionic and bosonic steady state current, respectively, which can serve as benchmarks for special cases of interacting many-particle dynamics. For fermionic currents, we identify a symmetry-induced enhancement mechanism in parameter regimes where the coupling between system and reservoirs is weak. This mechanism is broadly applicable provided the inter-particle interaction strength is small as compared to typical exchange interactions.

We establish a mathematically rigorous, general and quantitative framework to describe currents of non-(or weakly) interacting, indistinguishable particles driven far from equilibrium. We derive tight upper and lower bounds for the achievable fermionic and bosonic steady state current, respectively, which can serve as benchmarks for special cases of interacting many-particle dynamics. For fermionic currents, we identify a symmetry-induced enhancement mechanism in parameter regimes where the coupling between system and reservoirs is weak. This mechanism is broadly applicable provided the inter-particle interaction strength is small as compared to typical exchange interactions.

I. INTRODUCTION
Currents-the specific physical feature of non-equilibrium steady states of open systems subject to a potential gradient established by reservoirs-are a prominent topic of various branches of condensed matter physics. As the sizes of technological devices driven by currents reach mesoscopic scales, non-trivial quantum effects must be taken into account [1][2][3][4]. Yet, our theoretical understanding of currents in quantum systems is far from complete, e.g., many results are available for perfect lattices [5], but more realistic set-ups, with disorder and decoherence, still pose a panoply of open questions [6][7][8][9][10][11][12].
The past decade has seen a vivid debate on the relevance of quantum mechanics in biological systems and most notably in photosynthesis [13][14][15][16][17][18][19]. Since photosynthetic organisms are immersed in an environment of thermal photons, one may describe the situation via a constant influx of photons triggering an outflow of electrons [20]. The system, a large collection of intricately coupled chlorophyll molecules, is therefore constantly experiencing a flow of excitons which may be interpreted as a current. At present, the debate [21][22][23][24][25] on how such flow in the stationary state can be affected by quantum coherence on transient time scales remains widely open.
The aim of this contribution is to provide a rigorous mathematical treatment of currents in non-equilibrium quantum systems. To achieve this goal, we need a model which is analytically controllable. Therefore, we treat the coupling between the system and the particle reservoirs in a Markovian way, i.e. we ignore memory effects in the dynamics and use therefore a dynamical semi-group. Moreover, we focus on systems in which inter-particle interactions are sufficiently weak, such that the system can be described by effective model of free, i.e., non-interacting particles. For fermions, this implies that the shifts in energy levels associated to inter-particle interactions must be small compare to the energy-level spacings associated with the exchange interaction induced by the exclusion principle. In this scenario, we can derive bounds on the current, which are sufficiently tight to be saturated by properly designed systems. Our results thus also serve as a benchmark for studies of quantum transport in systems where interactions (or other non-linear effects) cannot be ignored. A violation of our analytically derived bounds is an unambiguous indicator of non-trivial interaction-induced effects, beyond mere many-particle interferences between indistinguishable particles.
To establish such a versatile theoretical approach which can handle the above diverse scenarios, and, in particular, also accounts for potential quantum statistical effects on transport, Section II of our present contribution provides a self-contained introduction to the mathematically rigorous framework of many-particle quantum currents. We will herein strongly rely on algebraic quantum statistical mechanics, a formalism which stems from mathematical physics. This algebraic approach to quantum mechanics of many-particle systems is indispensable to study infinitely large systems (as we do in Section V). Within this framework, we introduce a quantum version of the continuity equation, applicable to open system dynamics of the semi-group type [40,41]. We consider three contributions to the dynamics: a Hamiltonian part for the reversible particle dynamics, and two non-Hamiltonian parts which describe particle injection and extraction, respectively. For such non-equilibrium many-particle systems we derive several fundamental properties: In Sections III and VI we derive an upper bound for the particle current in the fermionic setting, and a lower bound for bosonic systems, respectively. In Section IV, we show that the fermionic upper bound can be saturated by appropriate design of the Hamiltonian part of the dynamics. The algebraic framework allows us to go beyond the standard Fock space formalism, which we illustrate in Section V, where we derive an upper bound for the current density in a ribbon, i.e., a 2D lattice system with shift-invariance in one direction, and a finite width in the other.
The strength of our contribution is that it makes no assumptions on the underlying single-particle Hamiltonian, and that it is applicable whenever the interaction between particles can be ignored to a good approximation. Hence, our approach does not only provide fundamental insight on the achievable currents in non-equilibrium quantum systems, but also opens novel perspectives for research in the fields mentioned above, where one may exploit the here identified design principle in a specific context.

II. MANY-FERMION SYSTEMS
We first provide an introduction to the algebraic formalism which describes many-fermion systems. The results presented in this section are well-known in the mathematical physics literature on quantum statistical mechanics [42][43][44][45]. In Sections III to VI, we apply this formalism to investigate the physics of currents in open quantum systems.

A. Fock space
It is common practice to describe many-fermion systems in terms of Fock space. This space is formally constructed using a single-particle Hilbert space H, also referred to as the mode space, as basic building block which provides all degrees of freedom of a single particle. As postulated by Pauli, identical particles are independent of labelling, a constraint which either leads to bosons or fermions. The wave functions of the latter species change sign under odd permutations of particles which is reflected in the fermionic n-particle Hilbert space (1) The anti-symmetrisation implies that the space H (n) is linearly generated by functions of the form Here S n denotes the permutation group of n objects, π a permutation, and sign(π) the signature of π. Note that these functions are generally not normalised and that they vanish whenever the single-particle wave functions are linearly dependent, as expected from fermions. Functions of the type (2) are often called Slater determinants. The fermionic Fock space Γ(H) constructed on H is built to accommodate any number of particles and therefore glues together all n-particle spaces: where C describes the vacuum component where no particles are present in the system. In the fermionic case and for a finite dimensional mode space H the direct sum breaks off: Fermionic Fock space can never harbour more particles than dim(H). Often it is assumed that Fock space is sufficient to describe general many-particle systems which is slightly inaccurate. Fock space can only accommodate finite numbers of particles. Both in the case where H is not finite dimensional, e.g., for infinitely extended systems (see Section V), or for a bosonic system (see Section VI), physics is much richer than Fock space. To study this larger realm of many-particle quantum physics, we must switch to a description in terms of observables and select a Hilbert space representation that matches the given physical situation.

B. Algebra of observables
The main tools at hand in Fock space are the creation and annihilation operators: a † (ϕ) and a(ϕ) respectively, where ϕ ∈ H. We work in a formalism of non-local creation and annihilation operators which have a straightforward interpretation: they create and annihilate a single-particle state ϕ. Their action is easily given on Slater determinants: where we have identified the n-particle vector The annihilation operator is the adjoint of the creation operator, its action on a Slater determinant is Indeed, as one may expect for fermions, sign bookkeeping is required. Fermionic creation and annihilation operators obey the canonical anti-commutation relations (CAR) These operators generate an algebra that forms the basic mathematical framework for the description of many-fermion systems with a given mode space H. The key idea of algebraic quantum physics is that the algebra of observables, rather than a Hilbert space, is the central mathematical object to describe large quantum systems. As a general algebraic framework and to contrast it with the Fock space representation above, we introduce abstract creation and annihilation operators c * (ψ) and c(ϕ) respectively, ϕ, ψ ∈ H. It must be emphasised that these objects are no longer linear operators on the Fock space, but merely generate a formal algebra determined by the basic relations The * is a formal operation which is the abstract version of the Hilbert space adjoint †. One then completes the algebra with respect to the unique C * -norm to obtain the C * -algebra A CAR of the CAR on H. [46] The completion is needed to apply general mathematical results and to describe dynamics in a controlled way. This framework is necessary to describe general many-particle systems with infinite-dimensional single-particle spaces; in these systems, we cannot describe all possible physics for all possible states (see Section II C) on the level of Fock spaces. In our present contribution, we strictly require this framework for the study of the quantum ribbon in Section V. In this formalism, observables are those objects O ∈ A CAR which are constructed using 1, c * (ψ) and c(ϕ) and which have the additional property that O = O * . In the context of many-particle systems, it is often useful to focus on polynomials of c * and c, in which contributions with specific particle numbers are related to definite orders. In this work, we focus solely on the simple class of single-particle observables corresponding to polynomials of order two.
A single-particle observable is essentially an embedding of an observable on the one-particle space into the manyparticle framework. In Fock space, one assigns a copy of the observable to each different particle in an additive way, e.g., the total energy of a system described by a single-particle Hamiltonian is the sum of the single-particle operators for each separate particle. The formal algebraic way to express this second quantisation is via the mapping Γ : B(H) → A CAR from the space of bounded operators on H to the algebra of observables, which acts as where we may select any orthonormal basis {η j } of H. In order to ensure that Γ(O) belongs to the algebra one has to impose the rather restrictive condition that O is a trace-class operator. It is not hard to check that different bases {η j } yield a same second quantised observable.
A specific example of interest in the discussion of particle currents for finite dimensional one-particle spaces H is the number operator N which literally counts the number of particles in the system. This operator is in essence of single-particle type, as it is given by Indeed, particle currents describe the in-and outflow of particles and therefore the behaviour of the observable N goes hand in hand with the behaviour of such currents. Not only does the algebraic formalism require a more abstract description of the observables in our theory, it also implies a more general structure for the quantum states which determine the statistics of measurement outcomes for these observables.

C. States
A quantum state is commonly associated either with a state vector ψ or with a density matrix ρ. Expected values of observables O are given by O = ψ, O ψ or O = Trρ O. This presupposes a specific Hilbert space representation of the physical system. The more general algebraic formalism starts with expectation functionals that allow for a probabilistic interpretation [47][48][49]. Thus a state is a linear functional ω : A CAR → C on the algebra of observables fulfilling the requirements These properties are respectively known as the normalisation and positivity conditions. A useful tool to describe states, and their perturbations, on a C*-algebra is the Gelfand-Naimark-Segal (GNS) construction [43,45,50,51]. This procedure associates a unique, but state-dependent, Hilbert space representation of the algebra to state ω. This representation returns the state as an expectation with respect to a state vector. Different states may, however, lead to inequivalent representations, which typically happens in the thermodynamic limit of many particle-systems. As an example one may consider Bardeen-Cooper-Schrieffer theory [52][53][54][55][56][57][58][59], where states with a finite particle density in the thermodynamic limit must be represented in a different Hilbert space than the Fock space which is constructed by exciting the physical vacuum (see Section II A). The GNS construction is a key result in algebraic quantum physics, which stresses that the properties of the system's state are essential prerequisites to study physical models.
In the present context, where only single-particle observables are considered, there are simple ways to describe the relevant expectation values. A notable fact is the existence of a linear operator Q ∈ B(H) for each ω, which serves as a (non-normalised) density matrix and is commonly interpreted as a covariance matrix: In the class of gauge-invariant quasi-free states, this operator Q suffices to fully determine the state. [60] In general, this is far from true and one can just say that Q characterises the single-particle statistics.
The fact that we are considering states on the CAR-algebra directly implies that 0 Q 1. The first inequality is necessary to fulfil positivity of the state, the second is a consequence of the fermionic behaviour and represents Pauli's exclusion principle. It follows [45] that for a general single-particle observable Γ(B) ∈ A CAR , with B a trace-class operator on H, This identity might not seem spectacular but it offers an enormous computational simplification. It is, therefore, one of the key ingredients in all the following sections of the present contribution.
If ω is a normal [61] gauge-invariant quasi-free state it can be shown that hence directly expressing the expected particle number in terms of Q which is now also a trace-class operator. This condition is also sufficient to guarantee normality.

D. Dynamics
In the spirit of the algebraic approach, it makes sense to consider the elements of the algebra A CAR as the dynamical objects in the theory, whereas the states remain unchanged at all times. This formally implies that we can consider a mapping Λ t1,t0 : A CAR → A CAR for an evolution from time t 0 up to t 1 . The first obvious requirements for a well-defined dynamics are These demands must be fulfilled for any choice of t 1 and t 0 . A more debatable [62][63][64] assumption on the dynamics is complete positivity which formally says that the system can be trivially embedded in a larger system without having to fear for loss of positivity. Such embeddings are also important to include internal degrees of freedom in the description. Complete positivity in other words guarantees that effective descriptions of only a subset of relevant degrees of freedom are possible. The formal mathematical phrasing requires an extension of the algebra by any matrix algebra M N to obtain A CAR ⊗ M N . We may now trivially extend Λ t1,t0 on A CAR to Λ t1,t0 ⊗ id N on A CAR ⊗ M N . When Λ t1,t0 ⊗ id N is a positive map for any N , the dynamics is said to be completely positive [65][66][67].
In addition to complete positivity, one may impose another demand which rarely holds exactly for a real physical system but often provides a very good approximation [68][69][70]: We impose a one-parameter semi-group structure on our dynamical map. The term "semi-group" implies divisibility of the map and hence the existence of a generator. Moreover, the generator is time-independent and thus the map is only governed by t = t 1 − t 0 . In other words, we can write the dynamics in terms of Λ t , and obtain that In general, we do not assume that the inverse exists, thus withholding the family of maps from being a full-blown group.
This type of dynamics is particularly useful due to powerful mathematical results. The results by Gorini, Kossakowski, Sudarshan [71] and Lindblad [67] are well-know, but only hold for algebras of observables which can be described by bounded operators on a Hilbert space. Nevertheless, Lindblad provided a more general recipe for completely positive, one-parameter semi-group dynamics on a C*-algebra A: He showed that any equation of motion of the type leads to a dynamical map with such properties. Hence, we may follow this prescription to engineer a dynamical system with the desired phenomenological properties. In other words, we do not microscopically derive a master equation but rather study one which has the correct phenomenology. In our present work, we follow and explore a model described by Davies [72]. From here onward, we assume that H is finite dimensional which is a considerable technical simplification. In section V, however, we will deal with translation-invariant systems and discuss how to cope with this more general situation. In particular, (18) allows us to write the generator of the dynamical semi-group in a form that nearly resembles the standard Lindblad form [72]: where α j , δ j ∈ H. To make sure that Ψ in (18) is a CP-map, one must impose [72] that θ in (19) is the * -automorphism determined by The Lindblad generators D a and D d describe the injection and extraction of particles into and from the system, respectively. With respect to the system degree of freedom, these terms mediate absorption and dissipation, thus the superscripts a and d. More specifically, fermions described by single-particle state vectors {α j } are injected into the system with positive rates {γ a j }, and particles which state vectors {δ j } are lost from the system with positive rates {γ d j }. Note that also temperature dependences can be accommodated within the positive rates {γ a/d j }. We consider systems of non-interacting particles, therefore, in accordance to (9), we must set h = Γ(H) with H = H † ∈ B(H).
We follow the model of [72], and many results of that paper are relevant for the present one. Specifically, we are interested in the dynamics of single-particle observables, given by x = Γ(B) with B ∈ B(H). Using (9), we insert Γ(B) into (19) followed by a straightforward computation [69,72,73] based on the anti-commutation relations (8), and we find that the relevant equation of motion for one-particle observables is given by That P, A, D ∈ B(H) directly follows from their definitions. Because the semi-group dynamics generated by (23) is Markovian, all rates must be positive, which in turn implies that A 0, D 0, and hence P 0. Moreover, for convenience in Section III B, we assume that A and D are strictly positive. For the bound on the current that will be derived later on, we can always consider the general case A 0 and D 0 using continuity.

E. Non-equilibrium steady states
Now that the equations of motion are determined, we observe that they can be solved exactly: We notice that, through its dependence on the absorption generator A, the second term is specifically related to the population of the system, via the particles that are pumped in. To infer the statistical distribution of measurement results associated with the observable Λ t (Γ(B)), we need to lift (24) to the level of states, by virtue of (14): An alternative perspective can be formulated by considering the object ω • Λ t ; because Λ t describes a dynamical map, it actually follows that, for any t > 0, ω • Λ t is a quantum state in its own right. In other words, we can treat the dynamics in the Schrödinger picture, by defining a family of states Intriguingly, equation (25) even provides us with an explicit expression for the Q(t) that appears in (14); by rewriting (25), we find where Q is the single-particle covariance matrix for the initial state ω.
Typically, at asymptotic times, pumped systems relax into a non-equilibrium steady state where finite currents are flowing. This limiting state has completely forgotten the initial conditions of the system. Put differently, generically each system observable converges to a multiple of the identity. The way to describe the asymptotic state, is by explicitly considering the limit t → ∞ in (24). To do so, note that since P > 0, generically, and therefore we find that Here NESS stands for non-equilibrium steady state. This,e.g., implies that the expected number of particles in the system converges to It is not hard to show [72] that the NESS state is the gauge-invariant quasi-free state determined by Q NESS .

III. CURRENTS
Non-equilibrium systems are typically characterised by the presence of currents, even when they reach a stationary state. In this section, we first discuss general properties of currents, determined by the "continuity equation" (33). We translate these results to a quantum mechanical setting to arrive at a sound definition of quantum particle currents. Finally, we extensively discuss a fundamental property of fermionic currents, which is one of our key results: The existence of a universal upper bound-irrespective of the specific potential encountered by the particle flow.

A. Particle Currents
We start by a formal definition of the particle current in the context of quantum master equations. The general structure in (19) presents us with the change of particles in the system over time. Because the number operator N = Γ(1), we insert B = 1 into (23) and obtain Note that the Hamiltonian contribution vanishes in the evaluation of (23) because [Γ(H), Γ(1)] = Γ [H, 1] = 0. This implies that the Hamiltonian, which is itself an observable of the form (9), conserves the total number of particles. We now study the particle current as a thermodynamic flux [40,41] and focus on its behaviour in the NESS (29). We note that, by definition of the steady state, the time derivative of the number operator is zero in the NESS, which yields, after combination of (31, 32) with (14), the balance equation We can now define the current of out flowing particles as where the absolute value is added because we focus on the magnitude of the current.

B. Bounding the current
Although expression (34) suggests that the current is independent of the Hamiltonian H ∈ B(H), it is implicitly present in Q NESS . Indeed, we can rewrite the current, using (29), to obtain In principle, this expression allows for a direct computation of the current, although this is generically intricate, e.g., when the operators in (35) do not commute. It is therefore instructive to derive a bound, to gain some general understanding of the parameter dependence of the current.
To do so, we first introduce the super-operator G can be split into a sum of two commuting terms, left multiplication by P − iH, and right multiplication by P + iH, respectively. Therefore, we may write where Generically, G is invertible and for positive definite P > 0 we can use the identity ds e −s(P −iH) X e −s(P +iH) .
Next, we compute from which it follows that We now introduce a symmetrised zero temperature Duhamel (or Bogoliubov) inner product [45,[74][75][76]: Here X and Y are general linear operators. Positivity of the scalar product follows from the invertibility of G, from G(X † ) = G(X) † , from (41) and from For an explicit evaluation of Schwarz' inequality we observe that from which we infer that Inserting these results in (44), we obtain and it then follows that which is the desired bound to the current. It is a universal one, since it does only depend on the reservoir coupling agents A and D, but not on the potential landscape which the fermions have to be transmitted through, defined by the system Hamiltonian H.
Since J max lacks a dependence on the single-particle Hamiltonian H, it is suggestive to inspect the tightness of the bound (50) for variable relative strength of unitary dynamics and reservoir couplings. For this purpose, we slightly rewrite (35) as where we introduced the parameter λ, to scale the relative strength of the Hamiltonian part of the dynamics as compared to particle loss and pump. λ → 0 completely cancels the coherent part of the dynamics while λ → ∞ makes the oscillatory Hamiltonian part much faster than the rates with which the systems couples to the reservoirs.
In the remainder of this section, we seek to numerically confirm the validity of (50) when λ in (51) is varied. To approach this problem, we consider many realisations of the system, each time choosing a random λ, random Hamiltonian H, and random channels A and D. For every realisation, the NESS current (51) is evaluated and compared to the upper bound (50). The results of this procedure are shown in Fig. 1, where specific choices for the random matrix ensembles were made: We consider a system of m modes, i.e., dim(H) = m. The Hamiltonian H is sampled from the Gaussian orthogonal ensemble (GOE) [77] which implies that it is a matrix whose entries are sampled from a normal distribution: ������� ������� The parameter v is related to the spectral radius (i.e., the largest eigenvalue in absolute value) and physically v/ √ m can be thought of as the typical (i.e., root mean squared) coupling strength between different modes. The matrices which describe the couplings between the system and the channels must be constructed so that they are always positive semi-definite. A standard method to generate random matrices fulfilling this constraint is to resort to the Wishart ensemble [78]. The latter ensemble is solely determined by the number of absorption (dissipation) channels, m A (m D ). For our numerical simulations, we set where W a and W d are m A ×m and m D ×m matrices respectively. They are generated by choosing random components according to The additional factor (m A + m D ) −1 in (53) is included to set the average eigenvalue of P (23) equal to 1. With this choice of ensembles, we can genuinely interpret λ (51) as the ratio of the frequencies of the coherent oscillations induced by H and the incoherent rates contained in P. Fig. 1 clearly shows that the bound (50) is valid for all realisations regardless of the magnitude λ. Nevertheless, we do observe that the bound is typically more accurate in the limit of dominantly coherent dynamics, characterised by λ 1. In this contribution, we will not attempt to understand the specific statistical properties which are obtained from the random matrix theory treatment. We do, however, note that in the limit λ 1 (or mathematically λ → ∞), the current is susceptible to changes in the Hamiltonian (recall eq. (51)) and that, therefore, a natural next step is to attempt to saturate the bound (50) in this regime. [79] IV. SYMMETRY ENHANCED CURRENT In this section, we investigate how an appropriate design of the system can generate a current close to J max (50). Because we are considering a designed system, it is reasonable to focus on the regime λ 1 where the coherent dynamics has a strong influence on the current. To get a maximal effect of the coherent contributions, we rigorously focus on the regime λ → ∞. This allows us to treat the problem using perturbation theory. In this limit, rapidly oscillating terms appear in (51) and by the Riemann-Lebesgue lemma [80] many contributions to J cancel.
The Hamiltonian can be represented in its spectral decomposition as Here the E k are the eigenvalues of H and R k are the orthogonal projectors on the corresponding eigenspaces of H. Using first order perturbation theory (where 1/λ is small), we compute If we want to saturate the bound on the current we have to design the R k in an appropriate way. Structuring Hamiltonians goes hand in hand with introducing symmetries. We therefore assume the existence of a unitary operator U , such that In order for such a symmetry to be useful, it must connect the couplings of the absorption channels A to those of the output channels D, leading to the requirement Given these additional structures, we can rewrite (56) as The fact that U and H commute, implies that U is block-diagonal with respect to the spectral decomposition of H This further implies which can in general not be cast in a more transparent form. However, in the case where the Hamiltonian H is non-degenerate (implying that, apart from (57), there are no unitary symmetries present in the system), we obtain that R k are rank-one operators. In this case we can express U as such that e iθ k are the eigenvalues of U . In turn this leads to By virtue of (56), where D is replaced by A, and (29), the right-hand side is exactly lim λ→∞ 2TrAQ NESS . Due to the symmetry (57), this implies that However, from the balance equation (33), we read that Both equations (64) and (65) can hold simultaneously only when TrAQ NESS = TrA/2, which implies that The second equality in (66) follows from (58). Inserting expression (58) into (50) directly yields which is exactly what we wanted to achieve. In words, we have shown that, in the absence of degeneracies in H, it suffices to find a unitary operator U which commutes with the Hamiltonian (i.e., a symmetry) and transforms D into A, in order to saturate the upper bound for the current in the limit λ → ∞. The most natural picture to associate with such a mathematical formulation, is that of a reflection symmetry. The limiting regime of λ can be seen as a rigorous way of demanding weak coupling, implying that the time scales of the system dynamics are much faster than those set by the rates with which the system couples to its reservoirs.
In realistic set-ups, this limit is never exactly achieved, therefore it is instructive to conduct numerical simulations to assess the deviations from the optimally achievable current, as a function of λ. This is done in Fig. 2: To make the simulation as general as possible, we start by sampling the unitary U introduced in (57) from the Haar measure [77]. [81] Because the matrix is unitary and random, we can always obtain a spectral decomposition We use this decomposition as the starting point for the construction of the Hamiltonian, which we define as So [H, U ] = 0 follows by construction. The choice of the uniform distribution for the eigenvalues E k is arbitrary, it simply serves to ensure that the typical level spacing-and hence the typical frequency of the coherent oscillations-is independent of the system size N . We sample A from the Wishart ensemble, see (53) and (54), but must take the constraint D = U † AU into account. When focusing on the regime of λ 1 in Fig. 2, we observe a similar trend as for Fig. 1. However, once we approach λ 1, we observe that, indeed, J ≈ J max for all realisations. The fact that the bound (50) can be saturated in the regime of dominantly coherent dynamics, λ 1, can be understood in a straightforward way: On the one hand, the rates with which the reservoirs couple to the system set time scales for particle exchange, which also governs the bound (50). On the other hand, however, coherent time scales, set by the Hamiltonian, determine how the particles explore the various modes inside the system. Therefore, if these coherent time scales are too slow, particles will linger in the modes where they entered the system where they block the path for additional particles due to Pauli's exclusion principle. Hence, the limit λ 1, guarantees fast redistribution of particles within the system, and, in addition, the design principle (58) and (57) guarantees a balance between input and output channels, such that particles can be extracted efficiently. In general, however, interference effects, incorporated in the fact that A, D, and H do not commute, make this naive picture more complicated. This is precisely why a mathematically rigorous treatment is important and non-trivial.
Furthermore, we stress that this discussion makes statements on the current J (51) relative to the bound J max (50). However, the maximal current J max itself depends on the absorption and dissipation channels (as governed by operators A and D). Therefore, when we rescale these parameters as A → γA and D → γD, it directly follows, from expression (50), that J max → γJ max . The results of Fig. 2 imply that for any such value of γ, the current J (51) can be optimised, for λ γ, by appropriately designing the system according to (58) and (57). However, because the value of the bound increases with γ, it is conceivable that a large value of γ (and thus large rates of particle exchange between system and reservoirs) can lead to large currents, even for slow coherent time scales, i.e., λ < γ. This hypothesis is, indeed, confirmed in Fig. 3. We can define the current J γ which results from rescaling A → γA and D → γD:  (51) is generated according to (69) with mode number m = 10 and random symmetry operator U from the Haar measure [77]. For each realisation, the absorption operator A (23) in (51) (52) and (69), for fully random and designed Hamiltonians, respectively. For the fully random systems, a single set of input and output channels, represented by A and D, respectively, is randomly chosen according to (53) and kept fixed. In the simulation of the designed systems we fix only A and generate D according to (58).
In Fig. 3 we see that typically the current increases when we increase the incoherent rates for particle exchange (by varying γ). However, we also observe that the bound is tighter in the regime of dominantly coherent dynamics as given by γ 1. In this regime, the designed systems give rise to optimal transport by saturating the bound, whereas we see fluctuations in the full random systems (see inset in Fig. 3). It must be noted that the double logarithmic scale of the plot masks these fluctuation.
Up to this point, we studied systems which contain a finite number of particles at all times. The strength of the C*-algebraic treatment and the formulation of the model in terms of a general CAR is the possibility to extend the setting to systems with an infinite number of degrees of freedom. In the following section, we consider systems that require a technical treatment based on current densities. We prove a generalisation of the bound (50) to a class of shift-invariant systems as commonly encountered in theoretical solid-state physics.

V. THE QUANTUM RIBBON
Above we focused on systems with a finite dimensional mode space, which excludes models with a translational invariance in some spatial directions. The latter situation requires to perform a thermodynamic limit, i.e., we first have to consider a finite subsystem, and subsequently perform a limiting procedure where the size of the system tends to infinity while the particle density remains finite [45]. We now consider such a model situation, with some inspiration from [82].
The specific system under consideration is a ribbon: A 2D system with translation invariance in one direction, and finite width in the orthogonal dimension. We assume that the system is accurately described by a tight-binding model and therefore the single-particle Hilbert space is given by a discrete lattice H := l 2 (Z) ⊗ C d , where d quantifies the finite width of the lattice. For k ∈ Z we denote by 1 {k} the sequence in l 2 (Z) with 1 at place k and 0 everywhere   (70) is generated according to (69) with mode number m = 10 and random symmetry operator U from the Haar measure [77]. A single absorption operator A (23) in (51) is drawn from a Wishart ensemble (53) with mA = 10, which is kept fixed for all realisations. The dissipation operator D is determined by the condition (58).
else. The mode space of our system can then also be seen as k∈Z 1 {k} ⊗ C d with the one-step shift along the ribbon given by 1 {k} ⊗ ψ → 1 {k+1} ⊗ ψ.

A. Shift-Invariance
Let us first focus on the space l 2 (Z). The Fourier transform F : l 2 (Z) → L 2 [0, 2π) can be defined through its action on the indicator functions 1 {k} ∈ l 2 (Z): A bounded operator operator A on l 2 (Z) is shift-invariant if and only if A = F −1 MâF . Here andâ ∈ L ∞ [0, 2π) . Therefore a bounded shift-invariant operator on l 2 (Z) corresponds to a multiplication operator on L 2 [0, 2π) by a bounded function on [0, 2π). Hermitian operators correspond hereby to real-valued functions and positive semi-definite operators to non-negative functions. This can straightforwardly be generalised to H = l 2 (Z) ⊗ C d : We say that a bounded operator X on l 2 (Z) ⊗ C d is translation-invariant along the ribbon iff X = (F −1 ⊗ 1)MX (F ⊗ 1), whereX : [0, 2π) → M d is a bounded matrix-valued function. If we now denote e k,l = 1 {k} ⊗ e l , with {e l } the standard basis in C d , we may write that e k ,l , X e k,l = 1 2π It also follows that, for two shift-invariant operators X and Y, It is useful to generalise (73) to the case whereX : [0, 2π) → M d is integrable. [83] In general, such a choice leads to an unbounded X.
To discuss currents, we are confronted with the problem that the global number operator N , which counts the number of fermions on the ribbon, is not an element in the CAR algebra over l 2 (Z) ⊗ C d . In fact, A does not contain any shift-invariant elements except for the multiples of 1. Shift-invariant elements are introduced by their local restrictions on finite subsets Λ ⊂ Z of the ribbon. To construct these local restriction, we define the appropriate projectors We can now consider the Λ-restriction Γ(P Λ XP Λ ) of "Γ(X)", which is a bona fide element of the algebra. Translationinvariance manifests itself by Γ(P Λ+1 XP Λ+1 ) being the one-step shift of Γ(P Λ XP Λ ). The global number operator corresponds to the choiceX(x) = 1 for x ∈ [0, 2π) and its restriction to Λ is just the number operator for the mode space l 2 (Λ) ⊗ C d , i.e., it counts the number of fermions on the compact domain defined by the restriction Λ.
Suppose that a shift-invariant 0 Q 1 determines the one-particle expectations (13) and that X defines a shift-invariant one-particle observable as above. Both Q and X are determined by matrix-valued functionsQ andX on [0, 2π) that satisfy the requirements that 0 Q 1 andX be real-valued and integrable. We can now consider the expectation of the density of Γ(X) in a specific state ω Q where we rewrite henceforth P n := P {1,...,n} where we introduce the "∼" to refer to densities in the system. Because of translation-invariance, there is no problem in fixing the leftmost site of the interval at 1. A small computation, similar to the type of computations used in proving Szegö's theorem [82], yieldsx

B. Currents in the Quantum Ribbon
The bound (50) on the fermionic current nicely fits with shift-invariance. For shift-invariant H, A, and D, both sides of the bound scale linearly with the length of the sub-interval on the ribbon that we consider. It then suffices to renormalise the inequality to obtain an analogous bound for densities.
The dynamics is a priori similar to the dynamics generated by (23), although we now specifically focus on the situation where H, A, and D are shift-invariant operators. We can simply repeat the arguments from Section II E above and obtain that in the long-time limit any state asymptotically converges to a shift-invariant state determined by the matrix-valued function The particle density ρ(ω) for a translation-invariant state determined byQ (77) is given bỹ Note that in our dynamical system we are typically dealing with a particle densityρ(ω) that changes over time.
We start by considering the evolution of the local number operator N n = Γ(P n ), which is described by (23): such that Note that, because P n is a projector, (14) yields ω Γ(P n ) = Tr P n QP n , and that (74) implies that commutators and anti-commutators of shift-invariant operators are again shift-invariant. Therefore, we may use (77) and evaluate where we already used that By definition of the non-equilibrium steady state and we therefore define the current density as For every x ∈ [0, 2π), we can apply (44) to find which can be rewritten and integrated to obtaiñ as a universal upper bound for the fermionic current across the quantum ribbon.

VI. BOSONIC SYSTEMS
Throughout the preceding parts of this contribution, we focused on systems of non-interacting fermions. Our methods are, however, also applicable to systems of non-interacting bosons. In this scenario, we must consider additional technical details related to the algebra of canonical commutation relations (CCR) [45,51,84]. One technical issue is that, for infinite dimensional mode spaces H, the bosonic algebra only allows us to define creation and annihilation operators in a representation dependent way. Another technical issue is that, even for finite dimensional H, states are not necessarily given by a density matrix on Fock space. Therefore, we here deliberately focus on systems with a finite number of particles, such that we remain in the Fock representation at all times.
The bosonic Fock space is defined (quite analogous to (1, 3)) as with The CCR can now we written in terms of non-local creation and annihilation operators b † (ϕ) and b(ϕ) which act on "Slater permanents" in a similar fashion as (4) and (5). These unbounded operators on Γ b (H) satisfy canonical commutation relations: In analogy to the fermionic case, we describe our dynamics in terms of the phenomenological master equation [85] d dt and, in analogy to (19)(20)(21) Again, {α i ∈ H} denote the single-particle state vectors in which particles are absorbed into the system, whereas {δ i ∈ H} denote the single-particle state vectors from which particles are dissipated out of the system. Because the creation and annihilation operators are unbounded, it remains to verify that this leads to a valid dynamical map, i.e., that it fulfils the conditions (16) and maps elements of the algebra onto other elements of the alegbra. To do so, we evaluate where the bosonic P is defined as with A and D as in (23). A fundamental difference between the fermionic and bosonic case is that the bosonic P in (93) is not necessarily a positive semi-definite operator on the single-particle space. We use (92) to evaluate the dynamics of a general single-particle observable where {η i } forms an orthonormal basis of the single-particle Hilbert space H. Straightforward integration of (91) leads to We now observe that the case where P is not positive semi-definite can lead to severe problems because it typically does not allow for the system to remain contained within the Fock representation at all times. This can be understood by assessing the time evolution of the particle number expectation value. We consider systems which are initially local with respect to the Fock representation, therefore the state is given by a density matrix ρ which acts on Γ b (H) and However, when P is not positive semi-definite, for generic ρ the asymptotic particle number is given by and thus diverges in the long-time limit. Physically this means that the system is unstable and never reaches a steady state. Therefore, we must impose that and therefore P 0, to ensure that the system remain confined to Fock space for all all times. This implies that systems of non-interacting bosons which absorb particles from an external reservoir require a sufficient (as quantified by (98)) amount of dissipation to ensure the existence of a well-defined NESS. [86] Having imposed condition (98), we find that the solutions to the bosonic and fermionic equations of motion, (95) and (24), respectively, are very similar, such that the same analysis as above can be repeated. The bosonic continuity equation is the same as the fermionic one when we write it in terms of P : In the NESS we find However, the definition of P has changed, so that the following balance equation between incoming and outflowing currents holds: This implies that we can still describe the current flowing through the system as J := 2Tr (DQ NESS ).
Remarkably, this definition of the current, together with (98), implies that we can next precisely follow the steps (36-48) of the proof for the fermionic bound. However, in (49) we did employ the explicit form of the continuity equation, and therefore this step differs from the present bosonic case. We now find that which implies and therefore The inequality (103) is remarkable because its derivation is largely analogous to that for the fermionic case, but it ultimately produces a very different phenomenology: There is no upper bound for bosonic currents in the NESS. However, bosonic currents are always stronger than a given quantity J min which is set by the channels. In systems where A comes close to D, while respecting (98), we see that the rate at which particles stream through the system can become arbitrarily large. Finally, we numerically scrutinize the lower bound (103). These results shown in Fig. 4 are obtained through evaluation of the exact expression for J: ds Tr e (−iλH−P )s Ae (iλH−P )s D , where P is given by (93). The parameter λ serves the same purpose as in (51) and Figs. 1 and 2, and the simulations are performed in a similar fashion as for Fig. 1: The value λ is chosen randomly in a way such that log 10 λ is uniformly distributed, whereas the Hamiltonians are sampled from the GOE (52). The choice of A and D is more subtle because of condition (98). To satisfy this constraint, we rather choose P and A from the Wishart ensemble (53) to subsequently determine D = P + A.
The results in Fig. 4 confirm the prediction by the lower bound (103) and show a drastically different behaviour compared to the fermionic case of Fig. 1. These results can be understood as a manifestation of quantum statistics. However, bosons do not disturb each other statistically when they start piling up (as happens when Tr(D − A) ≈ 0). Where the fermionic "repulsion" is often more important than the particle interaction, this is not the case for bosons. Hence, the assumption of absence of interactions for bosons is rather more stringent than it is for fermions. Therefore, one should be careful when interpreting these bosonic results when Tr(D − A) is small and particle densities become high.

VII. CONCLUSIONS
We described many-fermion and many-boson systems in which particles are incoherently pumped in and dissipated from the system, such that the total dynamics can be considered to be Markovian (memoryless). We prove that, in the absence of interactions between particles, the total particle current across the system exhibits universal properties in the stationary state: We could derive an upper bound (50) for fermionic currents, and, under some additional conditions which prevent the system from unlimited heating, a lower bound (103) for bosonic currents. Remarkably, both bounds are independent of the specific potential landscape the particles are transmitted through.
Numerically, we showed that, though counterintuitive, the bounds are typically sharp in the regime where the coherent dynamics' frequencies are high compared to the incoherent rates which determine the time scales of the reservoir coupling. This also led us to design Hamiltonians, as generators of the coherent dynamics, which saturate the bound in the limit where coherent dynamics is dominant. We proved that, in this limit, very general symmetry properties imposed onto the Hamiltonian suffice to achieve our goal. More specifically, we considered a unitary operator that commutes with the Hamiltonian and maps channels which connect the incoming reservoirs to the system onto channels which connect the system to the outgoing reservoirs. With these design principles, we can saturate our upper bound for fermionic currents. We note that the centro-symmetry [87][88][89][90], discussed in the context of optimal transport, is a special case of our present design principle. Hence, this work also improves the understanding of how such symmetries enhance quantum transport.
Our results offer a starting point for the investigation of several new questions, ranging from the relation of the here presented results to the Landauer formalism [2,91] to applications, e.g., in the quantum transport theory of disordered systems [92][93][94] or in the quantum statistics of non-equilibrium dynamical processes [95,96]. On a more fundamental level, the natural next steps are to investigate [97] how particle-interactions or other general sources of dephasing [9,10] impact the here derived universal bounds. In addition, it is a natural question to wonder what happens when the assumption of Markovian dynamics breaks down, e.g., it was recently shown [98], for a non-equilibrium spin-boson model, that the current is optimal for an intermediate coupling between system and reservoirs.