Giant rectification in strongly-interacting driven tilted systems

Correlated quantum systems feature a wide range of nontrivial effects emerging from interactions between their constituting particles. In nonequilibrium scenarios, these manifest in phenomena such as many-body insulating states and anomalous scaling laws of currents of conserved quantities, crucial for applications in quantum circuit technologies. In this work we propose a giant rectification scheme based on the asymmetric interplay between strong particle interactions and a tilted potential, each of which induces an insulating state on their own. While for reverse bias both cooperate and induce a strengthened insulator with an exponentially suppressed current, for forward bias they compete generating conduction resonances; this leads to a rectification coefficient of many orders of magnitude. We uncover the mechanism underlying these resonances as enhanced coherences between energy eigenstates occurring at avoided crossings in the system's bulk energy spectrum. Furthermore, we demonstrate the complexity of the many-body nonequilibrium conducting state through the emergence of enhanced density matrix impurity and operator space entanglement entropy close to the resonances. Our proposal paves the way for implementing a perfect diode in currently-available electronic and quantum simulation platforms.


I. INTRODUCTION
Due to their vast and exciting phenomenology, as well as the rapid advances on their control protocols, manybody quantum systems have become key ingredients for the development of novel technologies at the nanoscale.In particular, in nonequilibrium regimes such systems feature properties which make them highly appealing for applications in quantum circuits [1][2][3][4].This has been exemplified in several platforms where tunable transport of particles or heat can be induced and characterized, including setups based on electronic devices such as molecular electronics [5] and quantum dots [6,7], or in quantum simulators such as cold atoms [8,9] and superconducting qubits [10,11].In parallel to these seminal experimental advances, there have been numerous recent theoretical proposals to use quantum systems of interacting particles as different types of circuit elements that efficiently perform specific tasks.This includes autonomous quantum thermal machines [12][13][14], transistors [15,16], magnetoresistors [17], and a quantum analogue of the Wheatstone bridge [18].
Given their broad applicability, much attention in this field has been directed towards quantum diodes.For such devices, spatial asymmetries and non-linearities are engineered together to allow transport in one direction under a chemical potential or thermal bias, and to suppress it in the opposite direction when the bias is inverted.Experimentally, commercial electronic semiconductor diodes have demonstrated a ratio of forward and backward currents (or rectification coefficient) of ≈ 10 5 − 10 8 [19].
Here we put forward a rectification scheme which naturally emerges over a broad range of parameters in tilted interacting quantum lattices.These systems are the object of intense research as they have been argued to feature disorder-free (Stark) many-body localization (MBL) [37][38][39][40][41], quantum scars [42,43], and counter-intuitive phenomena such as time crystals [44] and transport opposite to an applied electric field [45].Moreover, they have been implemented experimentally in several quantum simulation platforms [46][47][48][49][50][51], making them attainable for state-of-the-art applications in quantum technologies.Neither the spatial asymmetry of the tilted onsite-potential nor the non-linearity of inter-particle interactions can on their own induce rectification.However, the ability to simultaneously engineer both ingredients in such platforms opens the possibility of its implementation.Here we reveal the underlying mechanism of this interaction-tilt interplay and show how the vastly differing transport properties [52] with the direction of bias give rise to giant rectification.Specifically, we find that it can exceed coefficients of 10 9 in small systems with no need of fine tuning and even with moderate interactions.We emphasize that in spite of the simplicity of our scheme, it features a rich phenomenology resulting in a genuine many-body rectification response, while also establishing a framework to build more complex nonequilibrium physics on top of it.
This article is organized as follows.In Sec.II we describe the model of our quantum diode, discuss how the transport properties are calculated and illustrate the emergence of large rectification in a small system.In Sec.III we demonstrate the insulating nature of the system for reverse transport.We explain the mechanism underlying current resonances for forward transport in Sec.IV.The combination of both results leading to giant rectification is discussed in Sec.V. Finally, in Sec.VI we present the conclusions and outlook of our work.Technical details of our calculations are described in the Appendices.

II. GIANT RECTIFICATION SCHEME
Our quantum diode is based on a boundary-driven configuration, which can be biased from right to left (reverse driving) or left to right (forward driving).For reverse driving, sketched in Fig. 1(a), both interactions and tilt each support on their own the formation of a high-energy current-blocking particle domain at the right boundary.When interactions and tilt coexist they thus cooperate to induce a high-energy state with an enhanced insulating nature.For forward driving, depicted in Fig. 1(b), both strong interactions and tilt favour on their own the formation of an insulating domain at the left boundary.For interactions alone this state remains high-energy, while for the tilt alone it is low-energy (see Fig. 1(c)).By being located at opposite ends of the energy spectrum, interactions and tilt now conflict with each other when they coexist.When both are of similar order, their competition breaks the domains and thus allows particle conduction.This forward-reverse transport asymmetry makes the system behave as a quantum diode.

A. Boundary-driven model
To unveil this phenomenon we consider a simple generic model of interacting particles: a one-dimensional lattice of spinless fermions in a tilted potential.The Hamiltonian for N sites is where c † j (c j ) creates (annihilates) a fermion on site j, n j = c † j c j is the corresponding particle number operator, J is the hopping amplitude (we take J = 1 to set the energy scale), ∆ is the nearest-neighbor density-density interaction, and µ is the diode's chemical potential.We take the tilt strength E > 0 so the potential always increases from left to right.For the noninteracting case, this model shows Wannier-Stark localization [53].Furthermore, it was recently proposed that for strong enough tilt the interacting model still features nonergodic behavior in the absence of disorder, an effect known as Stark MBL [37,38] and which is still the subject of debate [54].
To drive the system into a nonequilibrium state, we incorporate high-temperature reservoirs with differing chemical potentials at its boundaries.We assume a weak coupling to the system (Born approximation), memory-less reservoirs (Markov approximation), and that the bandwidths of the reservoirs are much larger than those of the system, which leads to frequency-independent system-reservoir interactions (wide-band limit) [55].Tracing out the reservoir degrees of freedom leads to a Lindblad master equation for the reduced density matrix ρ of the system [56], with the dissipative superoperators L k (ρ) defined by jump operators L k as The driving is induced by boundary operators , where Γ is the coupling strength (taken as Γ = 1) and the driving parameter −1 ≤ f ≤ 1 establishes the bias [4,57,58].Namely, when f > 0, fermions are mostly created on site 1 and mostly annihilated on site N giving forward driving, while inverting the sign of f gives reverse driving.The particle current corresponds to the expectation value of the operator which directly arises from the particle number continuity equation [4].We focus on the transport properties of the nonequilibrium steady state (NESS), where the current across the system is homogeneous, i.e.J j = J , as seen from the continuity equation.In addition, we set |f | = 1 to consider maximal forward/reverse driving, and take µ = −E(N + 1)/4 so the system features charge conjugation and parity (CP) symmetry.This choice will help make the rectification mechanism transparent.Nonetheless, under the assumed wide-band limit for the driving, the transport is independent of µ [59] (see Appendix C).

B. Calculation of transport properties in interacting systems
The transport properties of small lattices (N ≤ 8) were obtained exactly by diagonalizing the full Lindblad superoperator, written as a 4 N × 4 N matrix.The NESS Giant rectification in tilted strongly-interacting lattice.(a) For reverse transport, particles are mostly injected from a reservoir on the right boundary of the system, and mostly ejected to a reservoir on the left.For interactions or tilt alone, the resulting state is of high-energy, due to the adjacency of particles in the former case, or due to the particles filling the top of the potential in the latter.An enhanced insulator is induced when both coexist.(b) For forward transport, particles flow from the left to the right reservoir.The tilted noninteracting state is of low energy as particles now fill the bottom of the potential.Thus tilt now competes with interactions, which favors conduction.(c) Sketch of the eigenstructure of a tilted interacting system.i. High-energy insulating eigenstate preferentially populated by reverse bias.ii.High-energy insulating eigenstate preferentially populated by forward bias and dominant interactions.iii.Low-energy insulating eigenstate preferentially populated by forward bias and very large tilt.iv.The path from the interaction-to the tilt-dominated insulator, indicated by the dashed line, features avoided crossings at intermediate regimes where occupation probability is transferred between eigenstates, enhancing the transport.(d) Forward (F) and reverse (R) particle currents in a lattice of N = 4 and ∆ = 5 as a function of the tilt.Currents have been rescaled by the maximum current J0 obtained with the same driving for a noninteracting homogeneous system (see Eq. (A15)).(e) The corresponding rectification coefficient R.
is the eigenstate with zero eigenvalue, reshaped as the 2 N × 2 N density matrix ρ.For larger systems the calculation of the NESS is challenging.For forward transport, we describe ρ as a matrix product state [60], and obtain the particle current using the time evolving block decimation method [61,62] implemented with the Tensor Network Theory (TNT) library [63,64].Here, the evolution of an arbitrary state was simulated until the current becomes homogeneous, indicating convergence.However, for strong tilts and interactions a very slow approach to the NESS was often encountered making full convergence across the system impractical.In such cases we applied the following strategy.From the application of the continuity equation to the boundaries of the system [4], it is possible to relate the particle current to the population of the first site [65], so Thus, converging the population of the first site leads to a direct calculation of the current.This shortcut allowed us to accurately reproduce the NESS current of small lattices justifying its application to longer chains.
For N > 8 and reverse transport, where even the previous shortcut is insufficient due to exponentially slow convergence to the NESS [57,58], we calculated the current from an analytical ansatz for insulating domain states.This is given by (see Appendix B) with perturbative domain eigenstates |Ψ D n of each n particle sector.The ansatz captures remarkably well the particle density profile allowing the current to be calculated from the population of the first site.Since the considered system sizes can have very small currents J < 10 −15 , exact calculations from this ansatz were per- formed in MATLAB using the variable-precision arithmetic (vpa) function to take 32 digits of precision.

C. Rectification in small tilted systems
The simple nonequilibrium scenario introduced at the start of Sec.II is so effective that large rectification can be observed even for very small lattices.For instance, in Fig. 1(d) we illustrate the widely differing response of a system with N = 4 and ∆ = 5 when inverting the direction of the bias.For reverse driving of fermions down the tilted potential, the current decays rapidly and monotonically with E. On the other hand, for forward driving of fermions up the potential, the current features much larger values at well-defined resonances with E. As shown in Fig. 1(e), this leads to rectification coefficients of up to O(10 3 ), a value that can be vastly enhanced in longer chains (see Fig. 5).
These properties are intimately linked to the presence or absence of a particular distribution of fermions across the system, as depicted in Fig. 2(a).For reverse driving, the particles are pinned to the right boundary forming a domain half the size of the lattice, leaving the left half essentially unoccupied.The domain Pauli blocks more fermions from entering the system, resulting in an insulating state.On the other hand, for forward bias near the resonant values of E, the population profiles do not feature this constraint, so particles are allowed to move across the lattice.In the following Sections we explain in detail the emergence of these nonequilibrium effects.

III. REVERSE DRIVING INSULATING NESS
We consider first the insulating state.Domain formation is a common feature induced by strong interactions and tilt on their own, as depicted in Figs.2(b) where E = 0 and 2(c) where ∆ = 0.In both these scenarios there is no rectification; transport for forward and reverse drivings are identical up to a direction inversion, so we focus our attention to the latter.
For strongly-interacting homogeneous (E = 0) systems the insulating domain-like NESS is a well-known phenomenon [4,57,58,65], reproduced in Fig. 2(b) with exact diagonalization calculations for small lattices.The system features a flat profile in the noninteracting limit, typical of ballistic transport.As ∆ increases so does the slope of the profile, until reaching extreme population values at the boundaries for ∆ > ∼ 1. Increasing ∆ even further results in the creation of the particle domain.
In contrast, noninteracting (∆ = 0) boundary-driven tilted systems have received much less attention [66][67][68][69] and their manifestation of expected localization has remained unexplored.We perform such an analysis obtaining exactly the site populations n j and particle current J for large systems by reducing the Lindblad master equation (2) to a set of linear equations (see Appendix A).We use a rescaled tilt V = EN , to keep the on-site energy difference between boundaries approximately constant.This way the density profiles of different system sizes collapse for each V , evidencing the emergence of domains for small lattices and large tilt, as well as for long lattices and low tilt (i.e. when the total tilt is larger than the Bloch bandwidth, V > ∼ 4 [68]).In addition, for large systems the particle current decays exponential with N and V , as shown in Fig. 7 of Appendix A.
While the physical origin of the interaction-and tiltinduced insulating states is different, they share the same underlying mechanism that remains applicable for tilted interacting lattices.Specifically, they arise from the interplay between the gapped eigenstructure of the system and the maximal boundary driving.For small hopping J ≪ ∆ or E, the highest energy eigenstate |Ψ D n of each n particle sector (e.g.eigenstate ν = 5 of Fig. 3 is maximal, which breaks the insulating nature of the NESS.Also, when getting away from avoided crossings, the NESS is again mostly described by a single eigenstate, suppressing transport. a particle can be injected at site N or ejected at site 1) are exponentially suppressed (see Appendix B).Thus the eigenstates |Ψ D n become exponentially close to obeying making them dark states of the driving.The overall NESS density matrix ρ of the system is very well captured by an analytical ansatz built from a statistical mixture of these darks states introduced already in Eq. (6).From this, we show that ρ is strongly dominated by the contribution , corresponding to a state with occupied and empty domains of equal size (e.g.near unit probability for eigenstate ν = 5 for reverse driving in Fig. 3(b)).These results directly lead to the insulating NESSs of equal empty and occupied domains for systems with strong interactions [58], tilt or both, whose currents monotonically and exponentially decrease as either increase.

IV. FORWARD DRIVING RESONANCES
Next we focus on the forward driving setup with f = 1.Here, strong interactions and tilt on their own would favor a NESS largely dominated by perturbative corrections to the left domain |11 . . . 1 N 2 00 . . .0 N 2 .However, this configuration is the highest energy state of the n = N/2 particle number sector (degenerate with the right domain) for an interaction-dominated scenario, while it is the lowest energy state for a strongly tilted lattice.Thus, increasing ∆ raises the energy of this state, while increasing E lowers it.For instance, fixing ∆ and sweeping through E (as in Fig. 3(a)) closes the energy gap from the initial NESS dark state, eventually forcing it to cross through the bulk spectrum of the system (see sketch in Fig. 1(c)).As a result, current resonances are induced.This effect is cleanly illustrated by the halffilled even CP symmetric sector with the same parameters of Fig. 1(d) as shown in Fig. 3(a).Crucially, since all the eigenstates have the same symmetry, the nocrossing theorem [70] ensures that only avoided crossings occur.Their location accurately pinpoints the current resonances.
To understand the role of the avoided crossings we evaluate the contributions to the particle current in the energy eigenbasis.Since there are no coherences between different symmetry sectors, the current is given by where s sums over all the sectors, ρ s µs,νs = µ s |ρ|ν s is the coherence between eigenstates |µ s and |ν s of sector s, and J s νs,µs (q) = ν s |J q |µ s is the corresponding current matrix element evaluated on some site q.From the probabilities ρ s νs = ρ s νs,νs and coherences ρ s µs,νs , shown in Figs.3(b) and 3(c) respectively, the following picture emerges.For very low tilts and strong interactions, the NESS is solely dominated by the high-energy left domainlike state, evidenced by its large probability.This eigenstate has an energy ) to lowest perturbative order (see Appendix D), so as E increases it dives down towards the lower eigenenergies which have a weaker dependence on the tilt.As it approaches another eigenstate they mix and probabil- µ,ν and matrix elements J 2e ν,µ (q = 2), which gives a contribution from the half-filled even CP symmetry sector to the total current (dashed line).(c) Impurity (I) and OSEE (S).The second peak in both quantities, located between the first and second current resonances, emerges from the high population of eigenstates ν = 1, 3 (see Fig. 3(b)).In all panels, the dotted vertical lines indicate the location of the forward current resonances, corresponding to the peaks of the forward current depicted in Fig. 1(d) .
ity is transferred from the higher to the lower energy states.Thus, in the proximity of the avoided crossing the Hamiltonian no longer supports an eigenstate that is dark to the f = 1 driving.This allows large coherences to develop, enhancing the current.Further avoided crossings occur for even larger E, with lower eigenstates also getting populated.Eventually for very large E the whole spectrum is crossed and the lowest-energy eigenstate ν = 1 inherits the dominant role on the NESS.This leads again to an insulating left domain-like dark state, but now induced by the strong tilt.Notably, this picture is further evidenced by the behavior of the population profiles of the Hamiltonian eigenstates as the tilt is varied, discussed in detail in Appendix E. Overall, these results show that away from the resonances the eigenstate of highest probability at the NESS corresponds to a left domain, while close to the resonances the domains are suppressed and the eigenvectors involved in the corresponding avoided crossing exchange the form of their population profile.Some important remarks are in order regarding the transport resonance mechanism underlying our quantum diode.First, such mechanism is not restricted to half filling; different particle number sectors feature avoided crossings corresponding to forward current maxima, which might even be absent from the half-filled sector.Considering this, we have accounted for every resonance reported in this work.Second, in addition to being associated to a large coherence, a forward current resonance takes place at an avoided crossing if the corresponding current matrix element J s νs,µs (q) is significant.Otherwise, the contribution of that coherence to the current is small.This is already evidenced from the coherences shown in Fig. 3(c), and the current matrix elements of Fig. 4(a).The coherence between states µ, ν = 3, 4 has a large current matrix element at the tilt of the first resonance; the same holds for the coherence between states µ, ν = 2, 3 at the second resonance, and for that between states µ = 1, 2 at the third resonance.Because of this, the product ρ 2e µ,ν J 2e ν,µ (q = 2) is large, as depicted in Fig. 4(b), and a resonance is clearly seen.On the other hand, even though there is a large coherence between states µ, ν = 4, 5 for E ≪ 1, the corresponding current matrix element is strongly suppressed, leading to a barely-seen current maximum around E ≈ 0. In conclusion, not every large coherence, arising e.g. from an avoided crossing, will correspond to a forward transport resonance.
Finally, to emphasize the many-body complexity of the NESS at forward transport, we calculate the impurity I of the full NESS, and the operator space entanglement entropy (OSEE), [58,60,71], where λ α are the Schmidt coefficients of the half-size bipartition of the system.The former indicates how much the NESS deviates from a pure state.The latter has been well established as a measure of the complexity of quantum operators [72][73][74][75] (here a density matrix ρ of a dissipative many-body quantum system).Furthermore, being directly obtained from the singular value decomposition of the matrix product state representation of ρ, it quantifies the computational effort in capturing the nonequilibrium dynamics of the open system with tensor networks methods.Both quantities are plotted in Fig. 4(c), and feature maximal values at or very close to the current resonances.Intuitively we expect this to be the case: far from the resonances, where only one energy eigenstate is largely populated (so ρ is close to being pure), such a state approximately corresponds to a left product domain (low OSEE).Close to the resonances, the mixing of energy eigenstates naturally reduces the purity and also leads to the emergence of strong correlations, increasing the complexity of a matrix product state representation of ρ.These results reinforce that the observed enhanced conduction is a collective nonequilibrium phenomenon in an interacting many-body system.

V. GIANT RECTIFICATION
The combination of the enhanced insulating state of reverse driving in Sec.III with the forward resonant behavior of Sec.IV results in sizable rectification coefficients R = −J F /J R even for the small chain of the example (see Fig. 1(e)).Crucially, rectification becomes more substantial and more robustly achievable for larger system sizes N .The enhancement in the maximum R with N follows from the differing current scaling of reverse and forward drivings.Namely J R ∼ exp(−N ), characteristic of an insulator, while at a resonance we expect J F ∼ N −α with some exponent α > 0 characteristic of a conductor [2,4].As shown in Fig. 5(c), already for N = 8 with ∆ = 5 we find coefficients of up to O(10 7 ) compared to O(10 3 ) for N = 4.In addition, R is further enhanced by particle interactions.As depicted in Figs.5(a)-(c), while increasing ∆ diminishes both reverse and forward currents, the former decrease is much more prominent and thus their ratio features even higher resonances.Notably, here the current resonances also agree with maximal values of measures of complexity of many-body dynamics, as exemplified in Fig. 5(b) with the OSEE as a function of E for ∆ = 5.
While in small systems achieving the maximum R requires fine-tuning of E to locate a resonance, for larger N this constraint is rapidly removed due to the exponentially increasing number of eigenstates crossed within the bulk of the spectrum.As a result the rectification is broadened over the many-body bandwidth of the chain.Moreover, the favourable current scaling with N means that giant rectification can be achieved with a reduced interaction strength.To evidence these effects, using the insulating ansatz and tensor network simulations, we calculated reverse and forward currents for ∆ = 2 and systems of N > 8, shown in Figs.5(d) and 5(e) respectively.From these we find R ∼ O(10 9 ) over a range of E, as reported in Fig. 5(f).This establishes larger rectification coefficients than those previously obtained for different types of non-homogeneous lattices [32][33][34][35][36], with a more intuitive and easily implementable setup [47][48][49][50].
Finally, we note that the giant rectification effect discussed here is not an artefact of the driving scheme; it also emerges in more realistic boundary-driven setups.The Lindblad driving directly applied to the boundary sites of the system corresponds to fermionic baths with a flat and extremely wide spectral function at hightemperature [57].To move away from this idealized limit we explore in Appendix F a mesoscopic leads approach.
Here the infinite reservoirs at given temperatures and chemical potentials are replaced by leads consisting of a finite set of damped fermionic modes [12,68,76].We consider a simple but non-trivial driving scheme with leads comprising two modes that represent a double Lorentzian bath spectral function at high temperature and strong bias.As depicted in Fig. 11, even larger peaks of the rectification coefficient compared to those of Fig. 1 emerge, with both taking place at the same tilt values.Thus, we conclude that our giant-rectification quantum diode mechanism is implementable in more realistic nonequilibrium scenarios.

VI. CONCLUSIONS AND OUTLOOK
We have proposed a giant rectification scheme in correlated quantum systems, based on the natural asymmetry between forward and reverse transport in tilted interacting lattices.In our protocol, particles get locked into an insulating domain NESS when attempting to travel down the tilted potential (reverse bias).On the other hand, they propagate uphill resonantly as avoided crossings in the spectrum are approached (forward bias).Our calculations show rectification coefficients of up to O(10 9 ) for small chains, a number that can be pushed even further without fine tuning and at moderate interaction strengths by considering larger systems.
We emphasize that most of the reported rectification peaks are associated to sizeable values of the forward particle currents.This is demonstrated in Figs. 1 and 5 by rescaling J with the maximal current J 0 achievable under the same driving scheme for a noninteracting lattice with no tilt, which is a well-known size-independent ballistic current.Consequently at the resonances the system has sizeable conduction in the forward direction while it is insulating in reverse bias.Thus the observed giant rectification is not an artifact of dividing two very small currents associated with insulating states [31].
Our simple yet rich model serves as an ideal framework to unravel further nontrivial many-particle effects out of equilibrium.Namely, our results pave the way for developing different types of giant rectification protocols (e.g. of spin, charge or heat) in a variety of driving schemes, including finite temperature reservoirs, magnetic or thermal biases, and even excitation by light [77].For this purpose, recent approaches for accurately simulating more realistic reservoirs, even when they are strongly coupled to the system, might be exploited [12,68,76,[78][79][80][81].In addition, our work motivates further analysis on the nature of transport in the forward regime.This includes numerical simulations of larger driven systems, as well as analytical calculations of operators that overlap with the current based on the concept of dynamical l−bits [82,83].Furthermore, our proposal could be realized with existing electronic setups and quantum simulators.On the one hand, a tilted potential can be directly created on an electronic device coupled to metallic leads by an electric field [77].On the other, several experiments have implemented trapped-ion [46], cold-atom [47][48][49] and superconducting qubit [50] setups with a similar tilt and/or interactions to those used in the present work, and even in larger systems.Current quantum computing devices also constitute platforms where our scheme can be implemented [84,85].In addition, exploiting boundary-driven architectures in cold-atom platforms [9,24,86] or reservoir engineering for irreversible particle injection in superconducting circuits [87], similar versions of our diode could be incorporated in state-of-the-art nanoscale quantum circuits.In this Appendix we describe our approach to obtain the transport properties of non-interacting boundarydriven tilted systems.Even though we do not obtain the full solution for the density matrix, the relevant observables are calculated exactly up to numerical precision.For this, we propose the ansatz for the density matrix [59,88] where we define the l−site operators which correspond to generalized current and energy density operators respectively; in particular, B = −4J k .Even though this is a first-order expansion in f [59,88], it leads to an exact solution of the parameters a k /2.For this, the ansatz (A1) in inserted into the master equation ( 2) with the stationary condition dρ/dt = 0, which results in a closed set of linear equations for the parameters when grouping the coefficients in front of each operator.

Steady state equations
Importantly, due to the symmetric boundary driving several coefficients of the ansatz are zero or equal to others.Assuming N even, the symmetry relations are for the populations and 2 + 1 unknown coefficients, while for N odd we have (N −1) 2 2 + 1 unknown coefficients and N zero coefficients.This indicates that the same effort is required to solve systems of N = 2M and N = 2M + 1 sites.Considering these symmetries, we obtain the following set of equations after inserting the ansatz in the master equation and grouping the coefficients of each resulting operator.From the coefficients in front of 2n 1 − 1 we get This equation, consistent with Eq. ( 5), establishes a direct link between the boundary population and the particle current, allowing us to obtain one from the other.From the coefficients in front of B (2) This equation is modified for k = N 2 (i.e. when reaching the center of the system), since a N 2 = −a N 2 +1 and h (3) N 2 −1 .Thus here we have 2Ja N

+ 2Jh
(3) From the coefficients in front of For l ≥ 3 even, from the coefficients in front of B (l) and from the coefficients in front of where the δ lN results because for the last l the driving terms of both boundaries are equal and thus are summed.
For k = N −l 2 + 1, the corresponding equations are −2Jh N −l 2 +1 = 0. (A14) Finally, for l ≥ 3 odd, the same relations hold for k = 1, . . ., N −l+1 2 .This completes the set of N 2 2 + 1 linear equations for N even and of (N −1) 2 2 + 1 equations for N odd, which can be solved exactly for hundreds of sites with standard linear algebra packages and a moderate computational effort.Importantly, if E = 0, this problem reduces to that of a noninteracting ballistic conductor [59,88].Here the equations for coefficients in front of H for l ≥ 3 decouple from the rest, and form a set of homogeneous equations with the same number on unknowns; the solution of this system is the trivial one where the corresponding coefficients h k are zero.This leads to the size-independent particle current and a flat population profile in the bulk, as seen from Eq. (A7).In addition, if E = 0 but f = 0, the full system of equations is homogeneous, where only the trivial solution of all parameters being zero is possible.This corresponds to the NESS being equal to the identity (infinite temperature state), as expected.

Analytical results for small chains
Now we provide analytical results for the exact set of linear equations of ∆ = 0, focusing on the particle current and population profile for small lattices.We take first N = 4, which corresponds to a set of 9 equations and 9 unknowns, namely a 1 , a 2 , b (2) , h 1 and h (4) 1 .Even though this system can be easily solved, it leads to cumbersome expressions for the transport properties.Thus we consider two limits.First, for E ≪ 1 the particle current to lowest order in E is given by Thus the current decreases quadratically from the E = 0 result of Eq. (A15).A more interesting situation is that of very large tilt E ≫ 1.Here we get the current and the coefficients with a 3 = −a 2 and a 4 = −a 1 .This indicates that a large tilt leads to an insulating state.The current is suppressed as a power law of E, as well as the deviations of the populations from (1 + f )/2 and (1 − f )/2 on the left and right halves of the chain, respectively.Indeed, if f = 1, as in our numerical calculations, the populations n k = (a k + 1)/2 of the chain are (A19) This corresponds to the formation of almost fullyoccupied and empty domains on the left and right halves of the system respectively.If f < 1 the system is also insulating, with domains of population values closer to 1/2, as shown in Fig. 6 for larger systems.
In addition, we calculated the exact solution of the linear set of equations for N = 5, which has the same number of unknowns.In the limit of low tilt E ≪ 1, we get It has a very similar form to that of N = 4, but with larger deviation from the ballistic result.On the other hand, in the limit E ≫ 1, we get the current and f = 1 populations (A22) Thus in the large E limit the current for N = 5 is four times lower than that of N = 4, and the populations are closer to 1 on the left and to 0 on the right than those of N = 4.This means that even though the observables are of the same order in both sizes, the larger system favors more the formation of particle domains.Finally, it is important to note that these results already show that there is no rectification when ∆ = 0, as the current is ∝ f .Interactions are thus required in the present scheme for it to behave as a diode.

Transport in noninteracting tilted lattice
Using the exact solution described above, we calculate the transport properties for noninteracting tilted systems larger than those obtained analytically in the previous section.First we consider different values of f and E on lattices of N = 14, as show in Fig. 6.Similarly to the results for N = 4 of Eq. (A18), f just modifies the amplitude of the deviation of the populations from the value 1/2, leaving the overall form of the profile identical.The effect of E (and of V = EN , shown in Fig. 2(c)) is different, as it changes the form of the profiles.Indeed, as E increases these become more flat on both halves of the lattice, preventing the propagation of particles and thus driving the system to an insulating NESS.
Furthermore, considering Eq. (A6), the current is proportional to f , as in the analytical results of N = 4, 5. Given the latter, we analyze the current for the maximally-driven f = 1 case only, when varying V and N .First we evaluate J as a function of the system size for several values of V ; the results are shown in Fig. 7(a).The current decays monotonically both with V and N , as expected.Also, for large V > ∼ 4.5 an exponential trend with the system size is clear, indicating insulating behavior.For lower V the decay with N is much slower, or even not appreciable in the scale of the plot (as for V = 2); here we anticipate that the exponential fit will emerge on length scales larger than the system size.
We also study the behavior of the current as a function of V for fixed system sizes; the results are shown in Fig. 7(b).Importantly, for large-enough systems, we find an exponential decay of J with V from a value that decreases with N .These results evidence how rapidly an insulating state develops due to the tilt.

Appendix B: Ansatz for domain NESS
To unravel the structure of the insulating domain reverse-driving NESSs, consider first J = 0, where the Hamiltonian eigenstates are products in the Fock basis.The highest energy of the n particle sector corresponds to a single eigenstate with a domain on the right (only degenerate with the domain on the left when E = 0), denoted as |B n = |00 . . .0 N −n 11 . . . 1 n .Crucially, this state increases its separation in energy from the rest as E or ∆ increase, and is a dark configuration of the f = −1 driving [58].For low hopping J ≪ ∆ or E and right-toleft particle transport, the associated eigenstate |Ψ D n is perturbatively built from break-away configurations corresponding to the most inner particle moving to the left and the most inner hole to the right.These are denoted as |ψ n k p and |ψ n k h respectively for k hopping processes.Due to the energy gap ∆ + k 2 E between |B n and these configurations, the kth order correction to |B n in perturbation theory is ∼ |ψ n k h + |ψ n k p with amplitude Here we introduced the Pochhammer symbol y (k) = y(y+ 1)(y + 2) • • • (y + k − 1) = Γ(y + k)/Γ(y), and a factor f j different from 1 when the particle (hole) reaches site j = 1 (j = N ): In the noninteracting (∆ = 0) and homogeneous (E = 0) limits the amplitude reduces to with localization length scale ξ(Q) = 1/ ln(Q/J).Thus, in general the corrections to |B n are suppressed at least exponentially with k.In addition, since only the configurations with a particle on site 1 or a hole on site N are not Note that an identical description can be performed in the forward f = 1 driving scheme for the homogeneous interacting case (with degenerate left and right high-energy domains) and the noninteracting tilted system (where the left domain is the state of lowest energy).However, when interactions and tilt coexist this approach is no longer valid.From the dark states |Ψ D n we define the statistical mixture of Eq. ( 6), where the probabilities p n are determined from a detailed balance condition.Namely, if P n→m is the transition probability from particle number sector n to sector m due to the driving, the detailed balance condition is P n→n+1 +P n→n−1 = P n−1→n +P n+1→n , with probabilities The solution of these equations for an even system gives is fixed by the normalization condition in Eq. ( 6), and is also the highest of all the probabilities p n .Thus, the NESS of the system is largely dominated by the state with occupied and empty domains of equal size.
To provide a benchmark of the ansatz of Eq. ( 6), we reproduced the f = 1 analytical results of Eq. (A19).In particular, for any even N ≥ 4 and keeping only m = ±1 in Eq. (B5), the probabilities of the ansatz are, up to This leads to the general expression which agrees with the result of n 2 of Eq. (A19) for N = 4.We have also reproduced n 1 keeping up to O(E −4 ).

Accuracy of NESS insulating ansatz
Now we show that even for small noninteracting tilted systems, the aynsatz of Eq. ( 6) captures very well the observed ansulating behavior.For this we compare the values of the deviations of the populations from 1/2 (Fig. 8(a)) and the particle currents (Fig. 8(b)) obtained from the ansatz and from the exact solver for lattices of different sizes.Remarkably, even down to E = 1 the results from the ansatz are very close to the exact ones.In the presence of interactions, the comparison of both types of calculations has been given in Fig. 5(a) for reverse transport, again showing an excellent agreement.These results demonstrate the huge power of the ansatz to describe the nonequilibrium properties of a boundarydriven large-domain insulating NESS.
We emphasize that in spite of the similar description of the domain-like insulating NESSs induced by interactions and tilt on their own, there are some key differences.First, the former requires a minimun value ∆ = 1, as for lower ∆ the transport is ballistic [4,57,58,65]; the latter can be achieved with any tilt provided the chain is long enough.Second, as previously discussed, insulating domains emerge for E = 0 and strong interactions at f ≈ 1 only, while they can do it at any f for ∆ = 0 and long-or tilted-enough systems.

Appendix C: CP symmetry of tilted interacting model
Hamiltonian ( 1) with µ = −E(N + 1)/4 is invariant under simultaneous charge conjugation and centerreflection operations, e.g.U n j U † = 1 − n N −j+1 for the CP symmetry operator U .Also, since U 2 = I, its eigenvalues are ±1, corresponding to even and odd symmetry.However, only at half filling U commutes with the total particle number operator Ñ = j n j , since U Ñ U † = N − Ñ .Thus, only in this case we can consider particle number conservation and CP symmetry simultaneously.This property is used to build the Hamiltonian of the even CP half-filled symmetry sector whose energy eigenspectrum features the avoided crossings that coincide with the forward current resonances (see Fig. 3(a)).
Crucially, this does not mean that our giant rectification mechanism relies on the presence of CP symmetry.In fact, the particle transport remains identical when this symmetry is broken by taking µ = −E(N + 1)/4.This occurs because within the assumed wide-band limit, the energy gap between different particle number sectors (controlled by µ) is irrelevant.Only the energy differences between eigenstates within each sector, which are independent of µ, are important for our device.The CP symmetry just allows us to isolate even and odd eigenstates, transparently evidencing the connection between the avoided crossings and current resonances.In light of this discussion, also note that forward particle transport is identical (up to a direction change) to reverse transport with inverted tilt −E.

Implementation of CP symmetric Hamiltonian
To illustrate how the exact eigenstructure of the CP symmetric tilted Hamiltonian is obtained, we take N = 4.We initially consider the half-filling sector, where the particle number operator Ñ and CP symmetry operator U commute.First we identify each equivalent class (EC), corresponding to each set of basis states that are equivalent under CP symmetry, and their representative state (RS) [89].These are given in table I, along with the period of each EC.The basis elements of the even CP half-filled sector, which are common eigenstates of the number and CP symmetry operator, are thus The Hamiltonian of this symmetry sector, with the given order of basis elements, is Its spectrum is depicted in Fig. 3(a); due to the nocrossing theorem [70], it only shows avoided crossings.
On the other hand, the half-filled odd CP symmetric sector only has one basis element, namely with energy −∆/4.This completes the 6 states of the half-filling sector.
Beyond half filling the CP symmetry operator does not commute with Ñ .Thus we can build sectors of fixed number of particles or fixed symmetry, but not both.However, each case has redundant information.In the former scenario, the Hamiltonian of m particles is equal to that of m holes.In the latter, the Hamiltonians of even and odd CP symmetric sectors are also identical.To see this, note that the basis elements of each sector are the ± equal superpositions of the states of each EC listed in table II, namely The ± superposition states of 0 and 4 particles are disconnected from the other basis elements, and have energy 3∆/4.The states of 1 and 3 particles are connected with each other by the Hamiltonian, which has the matrix representation Similarly to the half-filled case, the spectrum of this Hamiltonian features avoided crossings only.This completes the 10 states of the sectors beyond half filling, and thus the 16 states of the N = 4 system.The above construction can be directly extended to larger systems.This is exemplified in Fig. 9 for the even CP symmetry sector of 10 sites, which corresponds to results presented in Fig. 5(d)-(f).

Appendix D: Perturbative calculation of eigenenergies
Given the many-body nature of our problem, it is not feasible to obtain exactly the eigenstructure of the Hamiltonian for large systems, even with the explicit use of the CP symmetry described in the previous subsection.Considering that our giant rectification mechanism manifests for strong interactions, we perform a perturbative calculation of the most important eigenstates in this limit.For this we consider first the case of zero hopping, where the eigenstates are product states in the Fock basis.Here, the energies for N sites and domains of n particles pinned to the left (L) and right (R) boundaries are Clearly both states are degenerate and have the highest energy of the n particle sector in the absence of tilt [58].For a finite tilt, E n,R is the highest eigenstate of the n particle sector, while E n,L is the second highest one for low E and the lowest one for very large E.
When turning on a very weak hopping J, these energies are corrected by second-order processes of the most inner particle jumping to the neighboring empty site and then jumping back.Assuming no degenerate states, the corrected energies are For N = 4, we identify Ẽ2,R = Ẽ0011 as the energy ν = 5 of the half-filled even CP eigenstructure of Fig. 3(a), which monotonically increases with E. Also, far from E = 2∆, Ẽ2,L = Ẽ1100 corresponds to the eigenenergy whose avoided crossings in the exact spectrum match the forward transport resonances, namely to ν = 4 for 0 A similar description holds for larger systems; this is seen in Fig. 9, where we plot Ẽ5,R and Ẽ5,L for N = 10 sites, indicated by green and blue dashed lines respectively.We also focus on the case of L = 4 to perturbatively calculate the other energies associated to the basis elements of the half-filled even CP symmetry sector.In the zero hopping limit, these are Turning on a weak hopping, and keeping up to O(J 2 ), the corrected energies are For tilts lower than those of the divergences (E = ∆, 2∆), these expressions work reasonably well, as seen in Fig. 3(a).In fact, the first crossing of Ẽ1100 and Ẽ1001,+ takes place at E ≈ 2.08, which is very close to the location of the first forward current resonance, E = 2.15.

Appendix E: Eigenstate population profiles
To further illustrate the mechanism underlying the forward transport resonances, we show in Fig. 10, for several values of the tilt, the population of the five eigenstates of the even CP symmetry half-filled sector discussed on Fig. 3.For a small tilt (Figs.10(a) and (b)), eigenstates ν = 4, 5 feature domain profiles pinned to the left and right boundaries respectively; the latter is always present and becomes more polarized as E increases, while the former is strongly modified by E. Namely, as the first resonance is approached (Fig. 10(c)), the left domain of ν = 4 is degraded and the population of ν = 3 starts developing a domain profile, which is fully established after crossing the resonance (Fig. 10(d)).A similar pattern of domain suppression and population exchange among the eigenstates involved in the corresponding avoided crossing takes place at the second (Figs.10(e) and (f)) and third (Figs.10(g) and (h)) resonances.Finally, at very large tilt it is ν = 1, the associated eigenstate of highest population, the one that features a left domain (Fig. 10(h)).Clearly, the absence of the left domain at the resonant values of E for the eigenstates ν of highest probability at the NESS (see Fig. 3(b)) is in line with our picture for forward transport.
Appendix F: Rectification beyond Lindblad driving: Mesoscopic leads approach Here we show that the rectification mechanism introduced in this work also emerges in more realistic nonequilibrium setups than that of Lindblad boundary driving.The latter corresponds to an idealized limit of infinite high-temperature reservoirs with flat and very wide spectral density [57].Now, we move away from it by employing the so-called mesoscopic leads description.This approach decomposes infinite fermionic environments at given temperatures and chemical potentials in terms of a finite number of lead modes, each of which is coupled to a residual bath described by Lindblad damping [12,68,76].In brief, we consider the tilted system of interacting spinless fermions, with the Hamiltonian H of Eq. ( 1), coupled at its boundaries to macroscopic reservoirs of noninteracting fermions at temperatures T α = 1/β α and chemical potentials µ α (h = k B = 1, α = L,R).This setup is depicted in Fig. 11(a), and described by the total Hamiltonian Here, H α is the Hamiltonian of reservoir α, and H Sα corresponds to its coupling to the tilted system; the latter is modeled as hopping processes from each reservoir mode to the closest site of the lattice.Each reservoir is characterized by its Fermi-Dirac distribution f α (ω) = e βα(ω−µα) + 1 −1 and its spectral density with ω α,m the energy of mode m of reservoir α, and λ α,m its coupling to the closest site of the lattice.Obtaining the dynamics and steady state of this boundary-driven system is extremely challenging.Interactions between its particles generally make the problem intractable from analytical methods; in addition, finite temperatures and strong coupling to the reservoirs can induce large non-Markovian effects.The mesoscopic leads approach has been recently shown to successfully overcome these complexities [12].Its key idea consists of reducing the macroscopic reservoirs to a finite collection of L lead modes, each one connected to its own Markovian bath, as depicted in Fig. 11(b).In this scenario, each mode contributes to the effective spectral density with a Lorentzian function: Here, ǫ α,k is the energy of mode k of lead α, κ α,k its coupling to the closest site of the lattice, and γ α,k its Markovian dissipation rate.By carefully selecting the number of lead modes L and their corresponding parameters [12,68], the original spectral density characterizing the system-reservoir couplings (Eq.(F2)) can be accurately approximated from Eq. (F3).
For simplicity, we calculate the NESS of an interact-ing tilted lattice of N = 4 sites, using a minimal but nontrivial realization of this updated driving setup with L = 2 modes per lead.This way, the NESS can be obtained with exact diagonalization.To mimic the conditions of Fig. Just as in the Lindblad driving scenario, an increasing tilt impedes reverse transport and leads to large resonances of the forward current.Furthermore, the latter agree with the avoided crossings of the eigenspectrum.Thus, the giant rectification mechanism reported in this work is not an artefact of the Lindblad driving, but generally emerges from the interplay between strong particle interactions, tilt and large bias.A full study of charge and thermal quantum diodes in this more realistic setup, considering larger leads and different driving conditions, will be performed in a future work.
FIG.1.Giant rectification in tilted strongly-interacting lattice.(a) For reverse transport, particles are mostly injected from a reservoir on the right boundary of the system, and mostly ejected to a reservoir on the left.For interactions or tilt alone, the resulting state is of high-energy, due to the adjacency of particles in the former case, or due to the particles filling the top of the potential in the latter.An enhanced insulator is induced when both coexist.(b) For forward transport, particles flow from the left to the right reservoir.The tilted noninteracting state is of low energy as particles now fill the bottom of the potential.Thus tilt now competes with interactions, which favors conduction.(c) Sketch of the eigenstructure of a tilted interacting system.i. High-energy insulating eigenstate preferentially populated by reverse bias.ii.High-energy insulating eigenstate preferentially populated by forward bias and dominant interactions.iii.Low-energy insulating eigenstate preferentially populated by forward bias and very large tilt.iv.The path from the interaction-to the tilt-dominated insulator, indicated by the dashed line, features avoided crossings at intermediate regimes where occupation probability is transferred between eigenstates, enhancing the transport.(d) Forward (F) and reverse (R) particle currents in a lattice of N = 4 and ∆ = 5 as a function of the tilt.Currents have been rescaled by the maximum current J0 obtained with the same driving for a noninteracting homogeneous system (see Eq. (A15)).(e) The corresponding rectification coefficient R.

FIG. 2 .
FIG. 2. Formation of particle domains.(a) Site population nj for N = 4, ∆ = 5 and the resonant tilt values, for forward and reverse driving.(b) Site population nj for a lattice of N = 8, E = 0 and different interactions ∆ for reverse driving.(c) Site population nx for ∆ = 0, different system sizes and rescaled tilt V , with rescaled position x = j/N ∈ [0, 1].

FIG. 3 .
FIG. 3. Transport resonance mechanism in interacting tilted lattice of N = 4 and ∆ = 5, illustrated in the even CP symmetry sector of n = 2 filling (sector s = 2e).In all panels, the dotted vertical lines indicate the location of the forward current resonances shown in Fig. 1(d).(a) Eigenenergies Eν as a function of the tilt E. For low tilts, the eigenstate ν = 4 corresponds to the corrected left domain state.The inset amplifies the first avoided crossing.The dashed lines correspond to second order perturbative calculations of the eigenenergies (see Appendix D).(b) Probabilities ρ 2eν of eigenstates ν on the NESS.The line of ν = 5, R corresponds to the highest energy eigenstate and reverse transport; all the others are for forward driving.Also note that the curve of ν = 1 is peaked at the first two resonances, even though it is not involved in an avoided crossing.This is because it corresponds to a perturbative correction of |1010 (see Appendix D), which is strongly coupled to the f = 1 driving.(c) Coherences ρ 2e µ,ν of pairs of eigenstates µ, ν for forward transport.At the avoided crossing of eigenstates µ and µ + 1, ρ 2e µ,µ+1

FIG. 4 .
FIG. 4. (a) Current matrix elements J 2e ν,µ (q = 2) for interacting tilted system of N = 4 and ∆ = 5.(b) Product of coherences ρ 2eµ,ν and matrix elements J 2e ν,µ (q = 2), which gives a contribution from the half-filled even CP symmetry sector to the total current (dashed line).(c) Impurity (I) and OSEE (S).The second peak in both quantities, located between the first and second current resonances, emerges from the high population of eigenstates ν = 1, 3 (see Fig.3(b)).In all panels, the dotted vertical lines indicate the location of the forward current resonances, corresponding to the peaks of the forward current depicted in Fig.1(d).

FIG. 5 .
FIG. 5. Transport and rectification in larger interacting tilted systems.(a) Reverse particle current JR as a function of tilt for N = 8 and different interaction strengths.The solid lines correspond to exact diagonalization results, and the darker dashed lines to currents obtained from our analytical Ansatz (see Appendix B).(b) Exact forward particle current JF for the same parameters.The dashed line represents the OSEE (divided by 10) for the system with ∆ = 5.(c) Corresponding rectification coefficient R. (d) Reverse particle current JR as a function of tilt for ∆ = 2 and different system sizes.Given that the analytical Anstatz slightly deviates from the exact results at low tilts E ≤ 0.5, as seen in panel (a) for N = 8, we obtain the reverse currents of N > 8 in that regime from matrix product state simulations.(e) Forward particle current JF for the same parameters.(f) Corresponding rectification coefficient R. The legends of system sizes spread over (d) and (e) apply to the three bottom panels.Also, the currents have been rescaled by the maximal current J0 obtained with the same driving for a noninteracting homogeneous system (see Eq. (A15)).
Appendix A: Transport in noninteracting boundary-driven tilted systems

FIG. 8 .
FIG. 8. Comparison of noninteracting NESS expectation values from exact calculations and insulating Ansatz.(a) Populations across a lattice of N = 10 and several values of E, for exact (solid lines) and perturbative (•) calculations.(b) Particle currents as a function of the system size; the same convention of colors and symbols of panel (a) is used.

FIG. 9 .
FIG.9.Eigenstructure of half-filled even CP symmetric sector for N = 10 and ∆ = 2.The solid lines correspond to the eigenenergies, and the dashed lines to the second-order perturbative corrections to the domains states pinned to the right (green) and left (blue) boundaries, namely Eqs.(D2).

∞
FIG. 11.Rectification with mesoscopic leads approach.(a) Scheme of tilted system coupled to infinite fermionic reservoirs at its boundaries.Inset: Fermi-Dirac distributions of left (L) and right (R) reservoirs for completely occupied and empty cases respectively.(b) Mesoscopic leads description, in which the infinite reservoirs are substituted by a finite collection of fermionic lead modes, each one connected to its own Markovian bath (blue squares).(c) Forward (F) and reverse (R) particle currents in a lattice of N = 4 and ∆ = 5 as a function of the tilt.We considered L = 2 lead modes per reservoir, with parameters TL,R = 10, µL = −µR = 100, ǫ L/R,k = ±0.5 and γ L/R,k = 1.Currents have been rescaled by the maximum current J0 obtained with the same driving (namely, equal parameters of the lead modes) for a noninteracting homogeneous system.The vertical dashed lines indicate the location of the maximal currents for the Lindblad driving scenario.(d) Corresponding rectification coefficient R. Inset: Effective spectral density (F3) of the discretized leads, with the contributions from each mode.
1(d), we consider ∆ = 5, high temperatures T L,R = 10, chemical potentials µ L = −µ R = 100 to induce large bias for forward transport (so the left (right) bath is almost completely full (empty), see inset of Fig. 11(a)), and µ L,R with opposite signs for reverse transport.For both reservoirs we consider lead mode energies ǫ L/R,k = ±0.5, and dissipation rates γ L/R,k = 1.The effective spectral density is the sum of two Lorentzians depicted in the inset of Fig. 11(d).The currents and rectification coefficient are shown in Figs.11(c) and (d) respectively.

TABLE I .
C1) Representative states for N = 4 half-filled system, states of the corresponding equivalent class and period.

TABLE II .
Representative states for N = 4 beyond half filling, states of the corresponding equivalent class and period.