Exact description of quantum stochastic models as quantum resistors

We study the transport properties of generic out-of-equilibrium quantum systems connected to fermionic reservoirs. We develop a new method, based on an expansion of the current in terms of the inverse system size and out of equilibrium formulations such as the Keldysh technique and the Meir-Wingreen formula. Our method allows a simple and compact derivation of the current for a large class of systems showing diffusive/ohmic behavior. In addition, we obtain exact solutions for a large class of quantum stochastic Hamiltonians (QSHs) with time and space dependent noise, using a self consistent Born diagrammatic method in the Keldysh representation. We show that these QSHs exhibit diffusive regimes which are encoded in the Keldysh component of the single particle Green's function. The exact solution for these QSHs models confirms the validity of our system size expansion ansatz, and its efficiency in capturing the transport properties. We consider in particular three fermionic models: i) a model with local dephasing ii) the quantum simple symmetric exclusion process model iii) a model with long-range stochastic hopping. For i) and ii) we compute the full temperature and dephasing dependence of the conductance of the system, both for two- and four-points measurements. Our solution gives access to the regime of finite temperature of the reservoirs which could not be obtained by previous approaches. For iii), we unveil a novel ballistic-to-diffusive transition governed by the range and the nature (quantum or classical) of the hopping. As a by-product, our approach equally describes the mean behavior of quantum systems under continuous measurement.


I. INTRODUCTION
Diffusion is the transport phenomenon most commonly encountered in nature. It implies that globally conserved quantities such as energy, charge, spin or mass spread uniformly all over the system according to Fick/Ohm's law where the diffusion constant D relates the current density J to a superimposed density gradient ∇n. Despite its ubiquity, understanding the emergence of classical diffusive phenomena from underlying quantum mechanical principles is highly non trivial. Early works based on field theory and perturbative methods [1,2] pointed out the possibility that interactions do not necessarily lead to diffusion at finite temperature, a question addressed then more rigorously by using the concepts of integrability [3]. These questions have then fueled many exciting discoveries in low-dimensional interacting systems [4]. A notable example is the ballistic-to-diffusive transition in quantum integrable XXZ spin chains [5][6][7][8][9][10], which also exhibit a superdiffusive point in the Kardar-Parisi-Zhang universality class [11][12][13][14][15]. These discoveries have motivated the generalized hydrodynamical descriptions of integrable systems [16,17], providing an elegant path to the question of diffusion at finite temperature [18], and paving the way to the description of diffusive phenomena based on perturbative approaches [19][20][21][22][23][24].
The out-of-equilibrium driving protocol illustrated in Fig. 1, where a system is coupled to external dissipative baths, has been crucial to unveil and characterize Figure 1. A stationary current J flows in a one-dimensional lattice when connected to left (L) and right (R) fermionic reservoirs, described by Fermi distributions f (ε) with different temperatures T or chemical potentials µ. The wiggly lines denote dissipative degrees of freedom acting on the system with rate γ. For a fixed difference of chemical potential δµ = µL − µR, dissipative terms are normally responsible for the Ohmic suppression of the current, J ∝ 1/N . .
In this paper, we develop a novel approach to characterize the bulk transport properties of quantum resistors which we show to be exact and systematic for a wide class of quantum stochastic Hamiltonians (QSHs). Our starting point is the Meir-Wingreen's formula [78,79] (MW), which expresses the current J of a system driven at its boundaries, see Fig. 1, in terms of single-particle Green's functions. We show that, for Ohmic systems, the MW formula supports an expansion of the current in terms of the inverse of the system size N . We illustrate how to perform practically this 1/N expansion, which reveals efficient to derive the diffusive current and the diffusion constant: we assume that, in the N → ∞ limit, diffusive lattices admit a simple description in terms of independently equilibrated sites and demonstrate that a wellchosen perturbation theory over this trivial state leads to the desired 1/N expansion.
We provide a comprehensive demonstration of the validity of our approach in the context of QSHs. Relying on diagrammatic methods and out-of-equilibrium field theory [80], we show that single-particle Green's functions of QSHs can be exactly and systematically derived relying on the self-consistent Born approximation (SCBA) -a generalization of previous results derived for a dephasing impurity in a thermal bath [50]. Equipped with this exact solution, and relying on MW formula, we explicitly derive the dissipative current flowing in the system and show that the Keldysh component of the single particle Green's function encodes the Ohmic suppression of the current. Then, we explicitly derive the asymptotically equilibrated state by "coarse-graining" of single particle Green's functions and validate our procedure to perform the 1/N expansion.
We illustrate the effectiveness and versatility of our approach for three different QSHs of current interest: i) the dephasing model [30-32, 81, 82]; ii) the quantum symmetric simple exclusion process (QSSEP) [35,[83][84][85][86][87] and iii) models with stochastic long range hopping [46,88]. The case studies (i) and (ii) illustrate the effectiveness of our approach, providing simple derivations of the current J and of the diffusion constant D, in alternative to approaches relying on matrix-product state [31,32,81], integrability [30] or other case-by-case solutions [33,35]. Additionally, we address previously unexplored regimes, by exactly solving the out-of-equilibrium problem with fermion reservoirs at arbitrary temperatures and chemical potentials. Our approach also allows to access twotimes correlators in the stationary state which were not described by previous studies. For case (iii), we show instead the ability of our approach to predict novel and non-trivial transport phenomena, namely a displacement of the ballistic-to-diffusive transition induced by coherent nearest-neighbor tunneling in one-dimensional chains. A by-product of our analysis is that all the results presented here apply also for system under continuous measurement, which are currently attracting a lot of interest in the context of measurement induced phase transition [42,44,46,88].
Our paper is structured as follows. Section II describes how the MW formula is a good starting point to build a systematic expansion of the current in terms of the inverse system size N . Section III presents QSH and shows the exactitude of SCBA for the computation of singleparticle self-energies. Section IV shows how our formalism allows to fully compute the transport properties of the dephasing model, the QSSEP and the long-range model. Section V is dedicated to our conclusions and the discussion of the future research perspectives opened by our work.

II. RESISTIVE SCALING IN FINITE-SIZE BOUNDARY DRIVEN SYSTEMS AND PERTURBATIVE APPROACH
In this section, we introduce generic tools aimed at studying diffusive transport in boundary-driven setups like those of Fig. 1. For these setups, the current is given by the MW formula [78]. In the simplified (yet rather general) situation, where the reservoirs have a constant density of states and the tunnel exchange of particles does not depend on energy, the MW formula reads (we assume e = = k B = 1): where f L(R) (ω) = [e (ω−µ L(R) )/T L(R) + 1] −1 are the Fermi distributions associated to the left and right reservoir with chemical potentials µ L(R) and temperatures T L(R) . G R/A/K are the retarded (R), advanced (A) and Keldysh (K) components of the single-particle Green's functions of the system. They are defined in time repre- , where the (curly)square brackets indicate (anti)commutation [89]. c j is the annihilation operator of a spinless fermion at site j. The Γ L(R) matrices describe system-reservoirs couplings.
Our aim is to establish a systematic procedure to compute diffusive current for large systems. The starting point will be the state of the system in the thermodynamic limit (N → ∞). By identifying in the MW formula (2) the terms leading to Fick's law (1), we motivate the simple structure of the problem for an infinite system size. In resistive systems, a fixed difference of density ∆n := n 1 − n N at the edges of the system enforces the 1/N suppression of the current (J ∝ ∇n ∝ ∆n/N ). It is thus natural to perform a perturbative 1/N expansion of the current on the N → ∞ state. We conjecture a possible perturbation scheme and show its validity in the context of QSHs.
Without loss of generality, we focus on discrete 1D lattice systems of size N [90]. In this case, the Γ L(R) matrices in Eq. (2) acquire a simple form in position space: [Γ L(R) ] j,k = Γδ j,1(N ) δ j,k . We also express the local densities in terms of Green's functions, namely The MW formula then acquires the more compact form: where we have introduced the local spectral densities and made use of the fact that dωA L(R) (ω) = 1.
The local spectral densities A L(R) (ω) exponentially converge in the thermodynamic limit N → ∞. This feature is generally expected and is illustrated in Fig. 8 for different classes of QSHs. This observation allows to establish that the 1/N scaling, proper to diffusive currents, must entirely arise from ∆n in (3). The possibility to ignore the size-dependence of the first term of (3) imposes strong constraints on the 1/N expansion of the difference of density ∆n in diffusive systems. If we write this expansion as one notices immediately that the leading term ∆G (∞) has to compensate the first one in (3), implying A sufficient but not necessary condition fulfilling this relation is obtained by imposing at each boundary: which will turn out to be satisfied for QSHs. These relations have a simple and interesting interpretation. In the infinite size limit, the flowing current is zero and thus the stationary value of the densities at the boundary can be computed by supposing that they fulfill a fluctuationdissipation relation or equivalently, that these sites are at equilibrium with the neighboring reservoirs.
Reinjecting (4) in the MW formula gives the current and as expected, we get the 1/N diffusive scaling. This relation tells us that the information about the diffusion On the left, spacial correlations in the infinite size limit are depicted. These decay exponentially as a function of the distance and are non-zero only within a finite length a. By coarse-graining the theory over this typical length, we obtain an effective theory (right) consisting of an ensemble of uncoupled sites with a finite self-energy at equilibrium.
constant is hidden in the 1/N correction to the density profile which can be in general a non trivial quantity to compute. However, we will see in the following that there is a shorter path to access it via the use of an infinite system size perturbation theory. The main idea of the 1/N perturbation is to find an effective simple theory that captures the relevant properties of the system in the N → ∞ limit. From there, transport quantities are computed perturbatively on top of this limit theory. To determine this effective theory, we conjecture that there is a typical length a beyond which two points of the systems can be considered to be statistically independent. Thus, by coarse-graining the theory over cells of size a, each cell becomes uncoupled and in local equilibrium, see Fig. 2.
The reasons motivating such factorization are twofold. First, the current is suppressed as 1/N in the large system size limit, so the infinite size theory should predict a null stationary current. Second, factorization of stationary correlations has actually been demonstrated for a certain number of diffusive toy models, most notably in the context of large deviations and macroscopic fluctuation theory [35,83,91,92]. For instance, it is known that the n th connected correlation functions of physical observables, such as density, generically behaves as N −(n−1) . Thus, it is natural to assume that for N → ∞, correlations must be exponentially decaying over a length a. We will show explicitly that in all of the examples studied, this factorization in the coarse-grained theory will turn out to be true and provide an analytic estimation for a in App.F.
We now put these assumptions on formal grounds. Let j andk be the spatial indices of the coarse-grained theory The relation between the different components R, A and K of the single particle Green's functions are assumed to describe uncoupled sites at equilibrium with a local self-energy Σj [80]. These conditions require then local fluctuation-dissipation relations of the form with retarded and advanced Green's functions which are diagonal in the coarse-grained space representation These relations fix entirely the stationary property of the system in the infinite size limit. The specification of the free parameters µj, Tj, ω 0 j and Σj have to be done accordingly to the model under consideration. We will see that they take a simple form for QSHs, namely the self-energy Σj is frequency independent and the µj, Tj ω limit can be taken taken in Eq. (9), as expected in the Markovian limit of the dissipative process [79].
To get the current, one needs to go one step further and understand which terms have to be expanded. The thermodynamic equilibrated theory does not exhibit transport, thus should be left invariant by the part of the Hamiltonian that commutes with the conserved quantity, for us the local particle density. It is then natural to conjecture that the perturbative term for the current is given by the dynamical part of the theory, that is, the part of the HamiltonianĤ dyn which does not commute with the local density. Thus, we conjecture that, at order 1/N , the current is given by : where the ∞ means the expectation value must be taken with respect to the infinite system size theory. This formula has the remarkable advantage that its computational complexity is very low since the coarse-grained theory is Gaussian. We remark that the 1/N expansion presented here is not a standard expansion in the hopping amplitude τ , since the latter has an exponentially large degenerate manifold of states at τ = 0.
In Sec. IV, we show explicitly how these ideas unfold for QSHs, by comparing computations done from the 1/N theory with the one obtained from the exact solution that we present in the following Section Sec. III. Understanding to which extent and under which conditions Eqs. (9,10) and (11) can be applied is one of the very challenging direction of study, in particular in the context of interacting quantum systems without bulk dissipative terms.

III. VALIDITY OF THE SELF-CONSISTENT BORN APPROXIMATION FOR QUANTUM STOCHASTIC HAMILTONIANS
In this section, we present a class of quantum stochastic models and associated Liouvillians (12), that describe either stochastic local dephasing or stochastic jumps of fermionic particles on a graph. The random processes are defined by a quantum Markov equation also known as a Lindblad equation. We will show explicitly two ways, exemplified by Eqs. (15) and (A1), to associate an underlying quantum stochastic model to such Lindblad equation, a process known as unraveling or dilatation [93][94][95]. Of particular interest for us is the description in terms of quantum stochastic Hamiltonians (QSHs) (15). It provides a way to resum exactly the perturbative series associated to the stochastic noise, which coincides with the self-consistent Born approximation (SCBA) for single particle Green's functions. This method was originally devised for the particular case of a single-site dephaser in Ref. [50] and we extend it here to more general situations. We will show in Section IV that, relying on SCBA, we can derive the diffusive transport properties of these models and show the validity of the assumptions underpinning the perturbative 1/N expansion presented in Sec. II.
Consider a graph made of discrete points, each corresponding to a site. To such graph we associate a Markovian process where spinless fermions on a given site can jump to any other site only if the target site is empty, see Fig. 3. We define γ ij ≥ 0 as the probability rate associated to the process of a fermion jumping from i to j and γ ji = γ ij the reverse process. The generator of such process is given by the Liouvillian, which acts on the density matrix ρ of the system: The total evolution of the density matrix ρ is in general given by where L 0 generates what we call the free evolution, in the sense that L 0 is quadratic in the fermion operators c i and the related spectrum and propagators can be efficiently computed with Wick's theorem [96,97]. Such Liouvillians can generally describe single-particle Hamiltonians or dissipative processes (coherent hopping, losses,. . . ). We will consider L(ρ) as a perturbation on top of this theory.
There exists a general procedure to see L(ρ) as the emergent averaged dynamics of an underlying microscopic stochastic, yet Hamiltonian, process. Lifting L(ρ) to this stochastic process is known as unraveling and there is not a unique way of doing so, see Fig. 3. The stochastic Hamiltonian can be treated as a perturbation  (12) which describes random quantum jumps between sites connected by an arrow. An arrow leaving and arriving at the same site represents a local dephasing. To a given Lindblad equation, we can associate multiple stochastic process (blue and green boxes), a process called unraveling (orange dashed lines). The Lindblad equation is recovered by averaging over the noisy degrees of freedom (full blue lines). We show that the unraveling in terms of quantum stochastic Hamiltonian (QSH) is particularly useful for the diagrammatic expansion of the theory.
in field-theory which requires the summation of an infinite series. Our strategy is to pick the relevant stochastic theory for which there exists a simple way to reorganize the summation and then take the average in order to get the mean evolution.
We now proceed to present the unraveled theory. Let dH t be the stochastic Hamiltonian increment, generating the evolution, which is defined by We work in the Itō prescription and consider stochastic Hamiltonians of the form W i,j t describes a complex noise and we adopt the convention that W ij * t = W j,i t . The corresponding Itō rules are summed up by Using the Itō rules to average over the noise degrees of freedom one recovers the Liouvillian (12). Finally, an other point we would like to emphasize concerns the connection to systems evolving under continuous measurements. Indeed, another way to unravel (12) is to see it as the average evolution with respect to the measurement outcomes of a system for which the vari- are continuously monitored and independently measured with rate γ i,j [93]. Although the physics is radically different at the level of a single realisation of the noise, on average it gives the same result than the prescription (15). Hence, all the results that will be presented for the mean behavior of our class of stochastic Hamiltonians also describe the mean behavior of systems subject to continuous measurements. The unraveling procedure corresponding to continuous measurements is described in detail in Appendix A.

A. Self-energy
We show now that the perturbation theory in the stochastic Hamiltonian (15) can be fully resummed, leading to exact results for single particle Green's functions. To perform this task, we rely on the Keldysh pathintegral formalism [80], which describes the dynamics of the system through its action S. The presence of dissipative effects can be naturally included in S using Lindblad formalism [98]. The action gives the Keldysh partition function Z = tr(ρ t ) where ψ = (ψ + , ψ − ) are Grassmann variables defined respectively on the positive and negative Keldysh time contours C ± . We follow the Larkin-Ovchinnikov's convention [99], in which the Keldysh action S 0 corresponding to the free-evolution L 0 is expressed in terms of the inverse Green's function G −1 namely All variables in the integral (18) are implicitly assumed to depend on a single frequency ω, which coincides with the assumption of stationary behavior, valid for our class of problems. The inverse Green's function G −1 is itself expressed in terms of the retarded, advanced and Keldysh green functions G R/A/K , defined in Section II: and whose diagrammatic representations in the time domain are given in Fig. 4. The causality structure of the Keldysh Green functions is enforced by the suppression of correlators ψ 2ψ1 = 0. This means that a retarded propagator can never become advanced, which pictorially translates into the fact that a solid line cannot switch to a dashed one.
The action corresponding to the Liouvillian term (12) reads [98] (20) which is a quartic action in the Grassmann fields. At the level of single particle Green's functions, the action S L is incorporated through the self energy Σ, defined as the sum of all one-particle irreducible diagrams. As in equilibrium field theory, the Dyson equation relates the full propagator to the bare propagator and the self energies Σ: To compute the diffusive current from MW formula, Σ must be know to any order; an a priori difficult task given the quartic nature of the action (20). Instead, rewriting the action at the stochastic level allows us to exactly derive the self-energy Σ and solve this problem. In the field-theory language, the unraveling procedure exemplified by Eq. (15) leads to the equivalent action where S sto is related to S L by the average E[] over the noise degrees of freedom: In formal terms, this transformation is reminiscent of a Hubbard-Stratonovich transformation where the action becomes quadratic in terms of the Grassmann variables. Note that the complexity encoded in Eq. (20) is preserved by the consequent introduction of the space and time dependent noise dW i,j t . However, the noise correlations imposed by Itō's rules (16) allow a dramatic simplification of the diagrammatic expansion in γ i,j of the Green functions within the stochastic formulation. Such simplified structure does not manifestly appear when working with the Lindbladian (averaged) formulation of the problem (20) (see Fig. 14 in Appendix B).
The resummation works as follows. In Fig. 5, we show the diagrammatic expansion of (21) up to second order in the stochastic noise γ i,j . The wiggly lines represent dW i,j t . Since we are interested in the mean behavior, we have to take the average over the noise degrees of freedom. This amounts to contract wiggly lines pair by pair. From the Itō rules (16), we see that upon contraction, a  The important consequence is that all the diagrams which present a crossing of the wiggly lines vanish because of the causal structure of the Keldysh's Green function, namely that G R (t, t ) is non zero only for t > t and conversely for G A . For a detailed proof of this statement, see Appendix B. In particular, the constraints of avoided wiggly lines establishes the validity of the self-consistent Born approximation (SCBA) for the self-energy of single particle Green's function and generalize the approach presented in Ref. [50]. SCBA allows a simple and compact derivation of all components as exemplified by the diagrammatic representation in Fig. 6. Namely, we have that in position space For the retarded and advanced components, this relation takes a particularly simple form since G

R(A)
j,k (t, t) = ∓ i 2 δ j,k in position space. Note that this simple expression is only valid when the two time indices are taken to be equal and comes entirely from the causal structure of the Green's functions in the Keldysh formalism. One way to see this is to evaluate the step function θ(t − t ) for the retarded and advanced Green's functions from the discrete version of the path integral presented in 9.2 of [80]. To get the Keldysh component G K , one has to solve the self-consistent Dyson equation : which is a problem whose complexity only scales polynomially with the number of degrees of freedom in the system (such as the system size N of the setup in Fig. 1). This solves the problem entirely at the level of singleparticle correlation functions. Remark that this applies to any model as long as the bare theory respects a Wick's theorem and its propagators are known. It allows a systematic study of quantum systems in the presence of external noisy degrees of freedom. This ability to calculate the Keldysh Green's function is crucial to give an exact description of out-ofequilibrium transport in dissipative systems, as we are going to show in the next section.

IV. APPLICATIONS
We now proceed to employ the self-consistent approach to showcase our 1/N expansion, presented in Sec. II, against a large class of QSHs that display diffusive transport.
The action describing the out-of-equilibrium setting represented in Fig. 1 has the form The first term in the action, S Bd , describes the exchange coupling with gapless non-interacting fermionic reservoirs of chemical potential µ L,R and temperature T L,R . The corresponding action, under the assumptions discussed in Section II, was derived for instance in Ref. [79]: where ψ a is a shorthand notation for (ψ 1 a , ψ 2 a ), L designates site 1 and R designates site N . The action S 0 is the quadratic action related to the intrinsic dynamics of the system, which can describe various situations from coherent dynamics to single-particle dissipative gains and losses [79]. In this paper, we will focus on one-dimensional nearest-neighbour coherent bulk hopping, which is described by the standard Hamiltonian, with τ the hopping amplitude. The corresponding action reads The free propagators are directly derived from the previous expressions of the action and read Notice that the reservoirs act, through the hybridization constant Γ, as natural regulators of the imaginary components of the non-interacting problem [80]. Finally S sto is the action corresponding to the QSH (22). As explained in the previous section, the demonstrated validity of SCBA for the Dyson equation (25) allows to derive exact expressions for the self-energies (24), and thus for the propagators of the full theory. Such solution allows to fully determine the transport properties of the system through MW formula (3). As shown in Section III, Equation (24) implies a particularly simple form for the advanced and retarded components of the self-energy: Importantly, in the geometry of Fig. 1, we can derive a compact and explicit expression of (25) for the diagonal terms G K (t, t) where we introduced the N -dimensional vectors and M is an N × N matrix with elements Notice that only G K carries information about the biased reservoirs, as can be seen from (35). The first term in (3) depends exclusively on spectral functions, which are readily derived from Eqs. (30) and (32), while Eq. (33) sets, through Eq. (4), the expression of the density differences at the edges ∆n.
Note that our analysis shows that the matrix M (36) is the key object encoding information about diffusion and it appears exclusively in the Keldysh component of the single-particle Green's function (33). A convenient way to understand this is to consider systems with singleparticle gains and losses that do not display Ohmic 1/N suppression of the current. It was shown in Ref. [79] that, while (32) remains valid in those systems, the matrix M Figure 7. Particular 1D discrete cases that will be of interest. Only the noise contribution is presented in this figure. In the dephasing model, all the sites are paired with themselves. For the QSSEP, the pairs are between nearest neighbours. In the long-range model, a given point is linked to all the rest of the lattice with a coupling decaying as power-law.
in (33) becomes 0 for these systems and the current saturates to a size-independent value. Thus, having a finitelifetime in the retarded and advanced Green's function is not sufficient to get diffusive transport. The imaginary contribution to the retarded/advanced self-energy, such as the one in (32), has the interpretation of a lifetime for the free single-particle excitations of the system, yet it is the Keldysh component of the self-energy that describes the consequences of dissipative scattering on the transport properties of the system. When M = 0, equation (36) gives us a linear profile for the density profile, which eventually leads to a 1/N diffusive contribution for the current as discussed in the II.
These considerations are those underpinning our general discussion about diffusive transport in Sec. II. We now turn to the case-by-case study of the specific QSHs depicted in Fig. 7. As said in the Introduction, we will focus on three one-dimensional models: the dephasing model, the quantum symmetric simple exclusion process (QSSEP) and models with stochastic long-range hopping. For the dephasing model, every single point on the lattice is coupled with itself by the noise. For the QSSEP, the noise couples each point with its neighbours. For the long-range, a given point is paired to all the rest of the lattice with a power-law decay as a function of the distance. These processes are illustrated for all three models in Fig. 7 and we will give more details about their physical motivations in the related sections.
Without loss of generality, in the oncoming analysis of the current J, we focus on a linear response regime in the chemical potential bias. We set an idential temperature for both reservoirs T L = T R = T and µ L → µ + δµ, µ R → µ − δµ. We expand Eq. (3) in δµ. One thus obtains, to linear order in δµ: where A(ω) is the edge spectral function, which coincides with A L/R (ω), because of the mirror symmetry of the class of QSHs that we will consider. The second term can be expressed in the form in which W is an N dimensional vector whose components are given by

A. Dephasing model
The dephasing model describes fermions hopping on a 1D lattice while subject to a random onsite dephasing coming from dissipative interactions with external degrees of freedom. In the language of Sec. III, this model corresponds to the case where all the points are paired with themselves, which results in substituting the rates in Eqs. (12) and (15) (see also Fig. 7). There are various limits in which this model can be derived. For instance, it can be thought as describing the effective dynamics of fermions interacting weakly with external bosonic degrees of freedom within the Born-Markov approximation [38]. In Refs. [31,32,81] it was shown, relying on matrix product operator techniques, that the dephasing model exhibits diffusive transport. Two-times correlators in the XXZ under dephasing was also studied in [52] and were shown to exhibit a complex relaxation scheme. For bosonic interacting systems, it was shown that the addition of an external dephasing could lead to anomalous transport [100,101]. Additionally, as discussed in Section III, the mean dynamics of this model coincides with the one where the occupation numbers of fermions on each site are independently and continuously monitored [45,102]. For this reason, the dephasing model has recently attracted a lot of interest as a prototypical model exhibiting a measurement rate-induced transition in the entanglement dynamics [43,44]. Finally, we note that in Ref. [30] a mapping between the dephasing model and the Fermi-Hubbard model was established. Although we will not discuss this mapping here, we stress that it implies that our method also provides the computation of exact quantities valid for equivalent systems governed by Hubbard Hamiltonians. The stochastic Hamiltonian for the dephasing model is readily obtained from the substitution (39), namely where B t denotes a real Brownian motion with Itō rule dB j t dB k t = δ j,k dt. The retarded and advanced selfenergies are obtained from Eq. (32) and read  while G R,A are obtained by inversion of Eq. (30) with inclusion of the self-energy (41). These functions are symmetric and given by, for i ≤ j [79,103]: where B 2 ) 2 − 4τ 2 /2. The related spectral functions at the system edges A(ω) = A 11 (ω) = A N N (ω) is represented in Fig. 8 for different system sizes N . It displays N peaks corresponding to the eigenspectrum of the system without dissipation. The width of the peaks is controlled non-trivially by the hybridization constant Γ and the bulk dissipation rate γ Dph . Plots for closely related quantities in the γ Dph → 0 limit can be found in Ref. [79]. In this nondissipative limit, the height of the peaks does not decay with the system size N . On the contrary, for γ Dph > 0, the peaks vanish in the N → ∞ limit, and the spectral function converges exponentially towards a smooth function A ∞ (ω) as shown in the inset of Fig. 8. One can analytically derive A ∞ (ω), as the retarded Green function (42) at the edges G R 1,1 = G R N,N converges to The exponential convergence of the edge spectral function is reproduced by all the other QSHs discussed below and verifies one of the preliminary assumptions exposed in Section II, identifying the density difference ∆n as the term entirely responsible for the 1/N suppression of the dissipative current in (3).
Our approach provides an efficient way to compute the second term in (37), through an explicit derivation of the matrix M : As we detail in Appendix D, the expressions (38), (42) and (44) allow the efficient derivation of the current (37) up to system sizes N 10 3−4 . As a consequence, we can systematically study the expected crossover from a ballistic-to-diffusive regime expected at length scales N * γ −1 Dph [31]. See also Appendix E for additional details.
Two main technical advances of our approach compared to previous studies [26,31,32,81,97,104,105] consist in its ability to naturally address reservoirs with finite temperatures T < ∞, accessing transport regimes left unexplored by previous studies and to access twotimes correlators in the stationary state. An important consequence of our analysis is that the rescaled conductance of the system, that we define as G = N J/δµ, has a non-trivial dependence on the temperature T and the dephasing rate γ Dph , namely In Fig. 9, we plot the coefficients (α, δ, η) across the parameter space (T, γ Dph ). From the plot, we identify three main diffusive transport regimes, R τ,T,γ , in which these coefficients are different. Note that the regions are not connected by sharp phase transitions but instead by crossovers, which appear sharp in logarithmic scale. Deep in the three regions, the rescaled conductance takes the approximate values In previous studies carried in the T → ∞ limit for the reservoirs, where they can be described as Lindblad injectors [79], the conductance G is assumed to be proportional to the bulk diffusion constant D [4,21]. The density profiles in the system (see App. E) clearly show that such interpretation cannot be extended to lower temperatures. The emergence of coherent effects between the system and its baths leads to finite-sized boundary effects, which do not allow the determination of the bulk diffusion constant through Eq. (46). To obtain the bulk diffusion constant we can use our approach to derive the density profiles inside the system and far away from its boundaries. We numerically verify Fick's law (1) in the bulk and find the diffusion constant to be which is double the conductance in the T γ Dph limit, as expected. At variance with the rescaled conductances (46), this quantity is not affected by any boundary effect and it is in agreement with previous analytical ansatzes, valid in the infinite temperature limit [31]. The independence of the diffusion constant (47) from the temperature at the boundaries is a consequence of the stochastic dephasing (40), which locally brings the system back to an infinite temperature equilibrium state regardless of boundary conditions. We thus see on this example that our approach allows to compute both the twoand four-points measurements of the resistance. Even for diffusive systems, the distinction between the two processes can be important.
To conclude our analysis of the transport in the dephasing model, we note that the different transport regimes in (46) explicitly depend on the stationary bias n 1 − n N , which suffers from boundary effects in some regions of the (T, γ Dph ) parameter space. We confirm with our exact numerical solution that this is indeed the case. This interesting bias dependence is beyond the scope of the present paper and left for future studies.

1/N expansion
Let us now show how the diffusion constant (47), that we obtained from our exact solution, can also be easily derived from the novel 1/N perturbative theory we introduced in Sec. II.
The first step is to fix the action of the infinite size theory S ∞ with the aid of the coarse graining procedure. We start by disposing the elements of G R/A/K i,j as a matrix and subdivide it in square cells of width a. We take the average over all the terms in the cell to obtain the effective Green function G R/A/K i,j , describing the correlations between theĩ andjth cell. This procedure is illustrated in Fig. 10-(right) for the retarded Green's function and increasing cell size (a = 1 corresponds to no coarse graining). As the cell size increases, G R/A/K i,j becomes a diagonal matrix with the off-diagonal terms vanishing as 1/a and exponentially suppressed with the distance |ĩ −j|.
This explicit calculation confirms the diagonal structure of G R/A/K and the reduction of the action to a sum of local commuting terms S ∞ = j Sj, where Sj is the action associated to thejth cell. To simplify the notations, we drop the tilde indices from now on and implicitly assume that the calculations are done in the effective coarse-grained theory. The diagonal terms of G R (ω = 0) are depicted in Fig. 10-(left) as function of frequency with G K shown in the inset. As a increases, the symmetry center of the functions changes to ω = −2τ converging to the black curves depicting Eqs. (9),(10). As mentioned before, the only free parameters that need to be fixed in the local theory are µ j , T j and Σ j (ω). For the dephasing model, we find that the self energy is simply given by iγ Dph /2. For a single site such an imaginary term was shown [79] to coincide with the effective action of a reservoir within the limit µ j , T j → ∞ while keeping the ratio µ j /T j fixed. Let n j be the local density at site j, . Interestingly, at leading order in 1/N , this relation turns out to be verified even at the microscopic level, i.e for a = 1. This tells us that the local equilibration condition of the infinite size theory is always true in our case. We furthermore suppose that in the coarse-grained theory, the expression of the retarded and advanced components will be given by a single-site two-level system, i.e we suppose the following expression for Sj: Where we absorbed the −2τ shift of frequencies in the integral. Expression (48) is valid in the bulk, independently from any value of µ, T at the boundaries. We check explicitly that the coarse-grained theory indeed converges a = 1 towards Sj as a is increased as shown in Fig. 10.
In the path integral formalism, the 1/N corrections to the current (11) is given by whereĴ is the current operator,Ĵ[ψ + , ψ + ] is the evaluation of this operator in the fermionic coherent basis on the + Keldysh contour, • ∞ := D[ψ ± ,ψ ± ]e iS∞ • and S dyn is the Keldysh action (29) associated to the contour integral ofĤ dyn defined in (11). Here we have explicitly thatĤ dyn = τ j c † j c j+1 + h.c. The current operator is in this case :Ĵ A straightforward calculation reported in Appendix C then leads to an explicit derivation of Fick's law: where ∇ is the discrete gradient ∇n j = n j+1 − n j . Equation (51), derived from the 1/N expansion, coincides with the exact result (47) in the whole parameter space. Such agreement validates the 1/N expansion as a systematic and efficient procedure to compute diffusion constants. From the computational point-of-view, note that the 1/N expansion did not resort to any numerical schemes and provided an exact expression of the diffusive constant, which could not be extracted explicitly from the Dyson equation (25).

B. QSSEP
In this section, we illustrate how our method can also be applied to the study of the quantum symmetric simple exclusion process (QSSEP) [35].
The QSSEP is a model of fermionic particles that hop on the lattice with random amplitudes which can be thought as the quantum generalization of classical exclusion processes [92]. Classical exclusion processes have attracted a widespread interest over the last decades as they constitute statistical models with simple rules but a rich behavior that is thought to be representative of generic properties of non-equilibrium transport. It has been particularly impactful in the formulation of the macroscopic fluctuation theory (MFT) [91], which aims at understanding in a generic, thermodynamic sense, macroscopic systems driven far from equilibrium. It is hoped that the QSSEP will play a similar role in a quantum version of MFT, which is for now largely unknown.
We are interested in a model of QSSEP plus the coherent jump Hamiltonian (28) that was first studied in Ref. [33]. The case of pure QSSEP can be retrieved in the limit τ → 0. As for the dephasing model discussed in Sec. IV A, we will see that the 1/N expansion formalism again offers a simple route to derive the diffusive current.
As pictured in Fig. 7, the QSSEP couples nearest neighbour sites. It is derived from Eqs. (12) and (15) by taking the prescription The associated QSH is (53) From Eq. (24), we get the advanced and retarded components of the self energies: (54) The retarded and advanced Green functions are given by inserting the bare propagators (30) and the self energy (54) into the Dyson equation (25) . These propagators can be directly derived from the ones of the dephasing model by making the substitutions γ Dph → γ QS and Γ → Γ−γ QS /2. As a consequence, all the considerations made for the spectral function and Fig. 8, in the dephasing model, equally apply to the QSSEP. This is not the case for the Keldysh component, where the M matrix has the different expression [106] Combining the above equation with (33) allows to obtain G K and allows to compute the current from (3), or its linearized version (37). For all values of the parameter space (T, γ QS ) the current follows the relation (see Fig. 11) which tells us that the diffusion constant is γQS 2 + 2τ 2 γQS in agreement with the result presented in [33]. For τ = 0, this generalizes the result from [35] which was restricted to boundaries with infinite temperature and chemical potential.

1/N expansion
The expression (56) for the current can also be obtained easily in the 1/N perturbative approach illustrated in Sec. II. The action in the infinite size limit is again of the form (48). From (54) we see that the expression of the self-energy is similar to the one of the dephasing model by simply replacing γ Dph by γ QS up to differences that tend to 0 in the infinite size limit. The current operator from site j to j + 1 in the bulk is given here bŷ The first part is easily evaluated to be −γ QS ∇n j /2 to first order in 1/N in the diffusive limit. For the second part, we simply need to redo the previous derivation by replacing γ Dph by γ QS . The term iτ (c † j+1 c j −c † j c j+1 ) then becomes − 2τ 2 γQS (n j+1 − n j ) which yields (56).

C. Long-range Hopping
Finally we turn to the model with long-range hopping from the noise (see Fig. 7). In this model each particle can jump to any unoccupied site with a probability rate that decays with the distance as a power law of exponent α. Power-laws appear naturally for instance in quantum simulation with Rydberg atoms [107][108][109] where they emerge because of the dipole-dipole interactions. Depending on the order of the interactions between atoms, different power laws can be reached. In the limit α → ∞, we get an "all-to-all" model, i.e there are random quantum jumps between all sites. These types of models have recently attracted interest as toy models to understand the interplay between quantum chaos and quantum information notably in the context of random unitary circuits [88,110].
For the long-range QSH we have (58) and the corresponding Hamiltonian is where N α = 2 N/2 k=1 k −α is a suitable normalization condition such that N ∞ = 2 and N 0 = N . The limiting cases of this model are the QSSEP and "all-to-all" model, respectively α = 0 and α → ∞.
For the long-range hopping the expression of the retarded(advanced) self-energy is As before, injecting the bare propagators (30), (31) and (60) in (25) yields G R(A) . As illustrated in Fig. 15 in Appendix C 3, this form of the self-energy is equivalent to the one derived for the dephasing model (41), with the only difference that the effective dephasing rate γ becomes site-dependent because of the presence of boundaries connected to reservoirs . We verified that the exponential convergence of the spectral function illustrated in Fig. 8, equally applies, as expected, for this model as well.
The M matrix is which combined to (33) yields G K .
In the absence of coherent hopping, there is a simple argument to conjecture a phase transition in the transport properties of the system at α = 3. If one considers the stochastic process (59) alone, its average has a simple interpretation as a classical Markov process, where the probability for a fermion at site 0 to jump to site j during a timestep ∆t, given that the target site j is empty, is p j := γLR Nα|j| α ∆t. For a single particle, this defines a random walk whose variance is given by v := j p j j 2 which is related to the diffusion constant via D = v/∆t. This diverges at least logarithmically for α ≤ 3. However, note that there is no simple reasoning to understand what happens if one were to study the model with the coherent hopping term as, a priori, a purely classical analysis does not hold anymore.
For the numerical computations, we fix γ LR = 1 and T = 1000 but the results are independent of the latter. In Fig. 12, we show the dependence of the linear response current with the system size for different values of α. When α is small, the current saturates in the N → ∞ limit, while for large values of α it decays as N −1 , as depicted in dashed gray line. This a signature of a ballistic-to-diffusive transition that occurs at a finite value of α.
To characterize this transition further, we look at the order parameter D −1 = − lim N →∞ ∇n/J. For diffusive systems, D −1 is the inverse of the diffusion constant and should be zero for ballistic systems. In App. E, we discuss the numerical fitting required to obtain D −1 from a finite-size scaling analysis. D −1 undergoes a second order phase transition at a critical power α c ≈ 2.87 (see the dark-blue dots in Fig. 13). When approaching the transition from the diffusive region, the diffusion constant diverges as D ∼ (α − α c ) 1.21 (see the gray dashed line in Fig. 13). It is quite remarkable and counterintuitive that setting τ = 0 pushes the diffusive regime to values of α < 3 instead of the opposite. A naive reasoning would suggest that the addition of a coherent hopping term would push the ballistic phase to values of α larger than the classical estimate (α = 3), as a finite τ would favor the coherent propagation of single particles across the system. We observe that the opposite is surprisingly true, and we leave the exploration of this effect to future investigations.

1/N expansion
For α > α c , the 1/N expansion is valid and we can compute D −1 in the limit of infinite temperature. The action in the infinite system size is again of the form (48) and the lifetime is fixed by (60).
Unlike the previous models, there is no simple analytic expression for the diffusion constant since its derivation depends on the system size. We provide a detailed derivation of the diffusive current in App. C. In Fig. 13, we depict the results of the 1/N expansion for various system sizes (full lines) and overlap them with the numeric solution of (37) (dots). Both methods agree up to machine precision which may be an indication that the 1/N perturbative approach is surprisingly exact even in the ballistic regime, α < α c .
As already highlighted above, the interplay between transport and coherence gives rise to a rich physics in the long-range hopping model, but understanding it in depth is beyond the goals of this paper and will be addressed in a subsequent work.

V. CONCLUSION
In this work, we provided a comprehensive analysis of the large system size properties of diffusive quantum systems driven out-of-equilibrium by boundary reservoirs. In particular, we showed that diffusive quantum systems can be described by an effective and simple equilibrated Gaussian theory, which allows a systematic way to compute their diffusive transport properties via an expansion in the inverse system size. We illustrated the correctness of our 1/N expansion by comparing to exact results we obtained, using a self-consistent Born method, for a large class of quantum stochastic Hamiltonians which show diffusive behavior. In particular, the self-consistent approach allowed us to explicitly derive the structure of the effective Gaussian theory, which consists of decoupled sites with a finite lifetime and where the effective equilibration and diffusivity is entirely encoded in the Keldysh component of local correlations.
As an illustration of the effectiveness of our approach, we computed the current in three models that have been of interest in the recent literature: the dephasing model, the QSSEP and a model with stochastic long-range hopping. For the dephasing model and the QSSEP, we illustrated the ability of our approach to extend the study of transport to situations with boundaries at finite temperatures and arbitrary chemical potentials. This allowed us to show how dissipative processes restore effective infinite temperature behavior in the bulk and explicitly derive the effective Gaussian theory via a coarse-graining procedure. For the long-range hopping model, our analysis unveiled that coherent hopping processes trigger diffusive behavior in regimes where transport would be ballistic in the exclusive presence of stochastic long-range hopping. This counter-intuitive phenomenon is a remarkable example of the non-trivial interplay between coherent and dissipative dynamics in open quantum systems, which could be efficiently addressed based on the self-consistent approach.
The validity of the self-consistent Born approximation for our class of stochastic Hamiltonians provides in principle the solution to the noisy version of any model whose bare action is Gaussian. Our proof is not limited by stationary behavior or by the one-dimensional geometry of the problems addressed in this paper, but can be extended to time-dependent and higher dimensional problems as well. This possibility opens interesting perspectives for the investigation of novel phenomena in a large class of problems. Extension of our approach could be devised to study quantum asymmetric exclusion processes [111][112][113], spin and heat transport, the dynamics after a quench, fluctuations on top and relaxation to stationary states and their extensions to ladder geometries or with non-trivial topological structure. These settings have been for the moment largely untractable, or were solved by case by case methods, for which we provided here an unified framework.
An important issue raised by our work consists in showing whether our description equally holds and provides technical advantage for studying the emergence of resistive behavior triggered by intrinsic many-body interactions with unitary dynamics, where the breaking of integrability leads to diffusive transport [1][2][3][4][19][20][21][22][23][24]. A priori, the arguments presented in Section II apply for any quantum systems which follows a local Fick's law and, as such, they have the potential for very broad applications. Additionally, it is commonly accepted that the phenomenology of diffusion is associated with integrability breaking and subsequent approach to thermal equilibrium [114][115][116][117][118]. Understanding if and how our approach can help make this link clearer is an exciting open question. In this respect, we also note that, because of the existing mapping between the Fermi-Hubbard and the dephasing model [30], the self-consistent Born approximation allows to compute exact quantities in the Fermi-Hubbard model. As far as we know, exact solutions for this model were only obtained in the framework of the Bethe Ansatz and it is thus interesting that a seemingly unrelated approach allows to obtain exact quantities as well. Whether a connection exists between the two approaches and whether the exact summation allows to compute quantities out of reach of the Bethe ansatz are interesting open questions.

ACKNOWLEDGMENTS
We thank L. Mazza for useful suggestions during the writing of the manuscript. This work has been supported by the Swiss National Science Foundation under Division II. J.F. and M.F. acknowledge support from the FNS/SNF Ambizione Grant PZ00P2_174038.
We also thank X. Turkeshi and M. Schiró for making us aware, in the final phase of the writing of this manuscript, of their work [105] before publication, where a study of the dephasing model from the point of view of Green's function has also been performed.

Appendix A: Unraveling to continuous measurement
In this appendix, we discuss the unraveling of Eq. (12) to a quantum stochastic differential equation describing a system under continuous monitoring. In the Itō prescription the stochastic equation of motion of a quantum system subject to continuous measurement of an observable O + O † at rate γ is given by [95] where L 0 describes the dynamics in absence of mea- ) and D O (ρ) = Oρ + ρO † − ρtr(Oρ + ρO † ). If we assume that at each link we have two independent measurement processes 1 and 2 with the same rate 2γ i,j and O 1,i,j := c † j c i and O 2,i,j := ic † j c i . The corresponding measured observables are , namely the so-called bond density and the current. It is straightforward to see that averaging out (A1), we get (12) again.

Dephasing model
For the dephasing model, the definition of the current in the bulk from site j to j + 1 is given bŷ The expectation value ofĴ j in the stationary state is given by where we used the Larkin rotation and removed the terms ψ 2ψ1 as they are always 0 for causality reasons.
Using the action associated to the coherent jump S τ (C3) we get, from (49), to leading order in 1 N : where ∞ means the average with respect to the bare action in the infinite size limit, where all the sites are uncorrelated. Using Wick's theorem and that ψ a jψ b j+1 ∞ = 0, the previous equation greatly simplifies : We can now use the bare action of individual sites (in presence of the dephasing noise): (C6) to obtain the explicit expression of the current from which we immediately read the diffusion constant D = 2τ 2 γ Dph .

QSSEP
For the QSSEP, the self-energy for an individual site is Σ j (ω) = γ QS − γQS 2 (δ j,1 + δ j,N ). The current in the bulk is given bŷ The first part of the current already scales like 1/N at order 0 in the S τ expansion. The second term is evaluated in the same fashion as for the dephasing model. This leads to and D = γQS 2 + 2τ 2 γQS .

Long-range hopping
For the long-range hopping model, the local current is defined from the local conservation equation of the particle number : Recall the expression of the self-energy at site j (60) : which is depicted in Fig. 15.
To get the current with the 1/N expansion, we take, as for the previous model, the 0 th order term in the first term in the expression of the current and the first order term in the second part. We obtain: For simplicity, we give in this paper only the expressions for the infinite temperature and chemical potential boundary conditions which amount to take Lindblad injecting and extracting terms (see [79]). The current at the boundaries is then given by: In the stationary state we have that ∀j ∈ [1, N ], J in j = J out j which leads to the following system of linear equation to solve in order to get the density profile : where n and v are N -dimensional vectors with elements n j and M is an N × N matrix such that and v j = −δ j,1 α L − δ j,N α R (C20) In this appendix, we present some important elements of the numerical implementation. The first step to compute any presented result is to stabilize and efficiently evaluate G R(A) (ω) at any ω. For the case of a uniform stochastic noise (e.g. free system, dephasing), a naive use of (42) would require evaluating the ratio of two polynomials of order O(N ), a notoriously difficult task for large N using floating point arithmetics. A possible solution would be to resort to arbitrary-precision arithmetic but this would entail a heavy speed cost.
We used for the results of the present paper the fact that G R(A) (ω) can be written as a ratio of polynomials and therefore, decomposed into a product of monomials To efficiently find the zeros and poles of G R [119], we note that the inverse of G R is a simple tridiagonal matrix with a generic form whose inverse is given by [120] T . Therefore, computing the poles and zeros of G R requires computing all the zeros of the sequences {φ i , θ i } L+1 i=0 , a task that can be done efficiently. If the matrix is invariant under a reflection along the anti-diagonal, it is enough to compute a single sequence instead, φ i = θ L+1−i . This is always the case in the models studied in the present paper. Since a i does not depend on ω, φ i is a polynomial of degree i with the initial conditions defined as φ 0 = 1 and φ 1 = ω + a 1 . One can efficiently find all the roots {z k } i k=1 of φ i using a Weierstrass-like recursive method [121,122], see Eqs. D3 and D4 for a second and fourth order scheme where W k is the Weierstrass weight. We chose these derivative-free schemes to avoid computing explicit derivatives that would slow down the computation. Choosing the correct initial condition is critical to the success of the scheme. To find the roots of φ i , we initialize the scheme with the roots of φ i−1 plus an extra root.  Figure 16.
Scaling of the current as a function of the system size in the dephasing model. From left to right: γ = 10 −3 , 10 −1 , 10 1 , 10 3 . As the dephasing increases, diffusion sets in at smaller system sizes. The vanishing dependence of J with the temperature indicates the crossover into the Rγ region (46).
We empirically found that the extra root should have a random position close to the middle root (after sorting by the real part) to guarantee the best convergence. This initial choice can still fail when some roots are located very far way from the others, which occurs for example for the model QSSEP. This happens when, at some step in the iteration, two roots coalesce and C (i) k diverges strongly. In order to stablize this divergence, we introduce a damping factor κ that suppresses large corrections z κ is a purely empirically value, which we typically take as κ = max(|b|). The role of κ is to slow down the algorithm and allow the coalescing roots to separate. Our root-searching algorithm has thus two parts: a quick search using a secondorder damped scheme, followed by a fourth-order damped scheme to precisely locate the roots. Once all the roots are recovered, we generate the new matrixT obtained from the estimates of the roots. We consider thatT is a good estimate only when max T −1 (0) ·T −1 (0) < 10 −10 . With the exception of the QSSEP, we find a typical value max T −1 (0) ·T −1 (0) ∼ 10 −13 for any system size.
Once the poles and zeros of G R(A) (ω) are computed, we proceed to compute G K using (33). To evaluate the M matrix, we resort to the residue theorem. If the poles of G R(A) are simple poles, the sum over residues can be computed in parallel only requiring the evaluation of the monomials {(ω − z k )}. We note that while each monomial (ω −z k ) is of order unity, a sequential multiplication can lead to overflown errors in the limit of large N . To avoid this problem, we multiply the monomials at random. If the algorithm fails to, within machine precision, separate two roots, the residue is computed from the contour integral instead.
The last step to compute G K and the current J, is to perform the frequency integral convoluted with cosh −2 ( ω−µ 2T ). This is done by evaluating the integral using a discrete integration scheme instead of residue theorem. Since the thermal dependence is only encoded in the cosh −2 ( ω−µ 2T ), discretizing the integral allows us deal with different (T, µ) values at no significant cost. We carefully verify that the mesh is fine enough to guarantee convergence of the integral at any (T, µ).

Appendix E: Finite-size scaling
In this section, we detail the finite-size scaling analysis necessary to plot Figs. 9,11 and 13.
The presence of a dephasing term is not enough to ensure that the system behaves diffusely at any system size.
Signatures of diffusive transport such as J ∼ 1/N , only emerge at a characteristic dephasing length, N * ∼ 1/γ. At short system sizes, or short time-scales, the system behaves as if it was ballistic. In Fig. 16 we highlight this ballistic-to-diffusive transition for different values of the dephasing and temperature in the baths. At small dephasing values, one cannot reliably extract the diffusion constant by fitting a a straight line to Fig. 16. Instead, In general, the diffusion constant decays with the inverse system size, which we exploit to extract the N → ∞ limit from a non-linear fit, dashed lines.
to extract the relevant information in the N → ∞, we use the fact that the diffusion constant has itself a 1/N scaling [21] when measured in the middle of the chain.
In the QSSEP and dephasing model, we use this result to perform non-linear fits to D as shown in Fig. 17. In this figure we plot the diffusion constant of the dephasing model as measured in the middle of the chain for increasing system sizes and different (T, µ) values. The dashed lines depict the non-linear fit of the function a + b/(N + c) with a, b, c fitting parameters. We find that most observables in these models exhibit 1/N corrections as discussed in [21]. The speed of convergence however depends on the point in the phase-space (T, µ), with region R τ (see Fig. 9) showing the slowest convergence. This is a consequence of the effects of the bath discussed in the main text. Deep in the τ -dominated regime, we observe the breaking of Fick's law near the edges as shown in Fig. 18.
Since this effect only occurs in a finite portion of the system close to the edges, the convergence is only slowed down. We thus evaluate D in the middle of the chain to mitigate its effects and get a better accuracy.
For the long-range model, one needs a different approach to obtain the N → ∞ limit correctly, especially when close to the ballistic-diffusive transition described in the main text. A tentative form for the finite size extrapolation is provided by the solution of the diffusion equation for single particle under a random walk with long-range hopping discussed in Sec. IV C which gives a diffusion constant D = H , correctly captures the finite-size dependence of D −1 for all α values. The fitting parameters a, b, c respectively describe the amplitude, critical exponent and possible finite-size corrections. In Fig. 19, we depict D −1 against the result of the fit, respectively dots and dashed lines. The best fitting parameters are plotted in the inset. The quality of the fit allows us to conjecture that, at the transition point, the diffusion constant diverges logarithmically D LR (α = α c ) ∼ H In this section, we analytically estimate the coarsegrain length a from the correlation length of the dephasing model. Due to Eq.(9), it is enough to estimate a from a single Green function, in this case the retarded component. The starting point is the analytic expression of the elements G R i,j in the bulk of the chain. For large systems, the boundaries become irrelevant and the good basis of the problem is the momenta basis. In k-space, the self-energy takes a diagonal form  For the QSSEP, there are cross-diagonal terms in momentum that vanish as 1/L and can be safely ignored. Since both self-energy and Hamiltonian are diagonal in the momenta basis, one has where k = 2τ cos(k) is the eigenenergy of the bulk Hamiltonian. To find the retarded function in position space, we take the Fourier transform with the continuum limit for k G R r,r = dk 2π e −ik(r−r ) ω − 2τ cos k + iγ/2 .
The integral can be solved using the residue theorem and, after some lengthy yet simple manipulations, we find a compact formula G R r,r = i |r−r |−1 2τ cos y e iy|r−r | , where y = arcsin ω+iγ/2 2τ is a complex variable with Im(y(ω)) > 0. Therefore, in the dephasing model an estimate for the correlation length is given by ξ = 1 min Im arcsin ω+iγ/2 In the limit of small dephasing γ, we have ξ = 4τ /γ which serves as an estimate for the coarse-grain length a ∼ τ /γ. As expected, a should be of the order of the dephasing length N * ∼ 1/γ.