Dissipation-Induced Mobility and Coherence in Frustrated Lattices

In quantum lattice systems with geometric frustration, particles cannot move coherently due to destructive interference between tunnelling processes. Here we show that purely local, Markovian dissipation can induce mobility and long-range first-order coherence in frustrated lattice systems that is entirely generated by incoherent processes. Interactions reduce the coherences and mobility but do not destroy them. These effects are observable in experimental implementations of driven-dissipative lattices with a flat band and non-uniform dissipation.


I. INTRODUCTION
The wave functions of a perfect crystal are described by Bloch states which are delocalised over the entire crystal. Transport through these states is possible as particles put into the system by a source can be extracted at a remote location by a drain. However, interference can strongly affect the crystal's transport properties. For example, disorder or impurities in the crystal create multiple scattering centers, leading to Anderson localisation [1][2][3]. Destructive interference localises the wave functions which can no longer transport particles through the crystal. Recently, the joint effects of interactions and interference resulting in many-body localization have attracted enormous interest [4][5][6][7][8]. Also in the context of biological systems, the effect of dephasing on interference in transport processes has recently been studied for light harvesting complexes [9][10][11][12].
Localisation can also exist in crystals without disorder where geometric constraints on the tunnelling rates lead to destructive interference. Frustration quenches the kinetic energy of the Bloch states resulting in a flat band with an infinite effective mass. This macroscopic degeneracy allows localised wave functions to be constructed which are insulating stationary states of the system. Synthetic flat band crystals have recently been created in photonic lattices [13][14][15][16], polaritons in etched semiconductor heterostructures [17,18], ultracold atomic gases in optical lattices [19] and surface plasmons [20,21]. Moreover, concepts for the realisation of such lattices have been put forward for superconducting resonators [22].
In several recently considered samples the explored particles of the system are photons [23]. Due to unavoidable experimental imperfections, these are usually prone to loss. In order to repopulate the crystal, a constant coherent or incoherent drive can be applied. Therefore, these driven-dissipative frustrated lattices provide a fascinating new platform with which to study localisation, coherence and mobility phenomena.
The first works that have recently started to explore geometrically frustrated lattices in such regimes investigated to what extent frustration phenomena of closed * e.owen@hw.ac.uk systems, such as the enhancement of interaction effects, survive in the presence of driving and dissipation [24,25]. Whereas the effect of frustration and interactions has been extensively investigated [26][27][28][29][30][31][32], the interplay of frustration and local dissipation is much less explored.
Here we show that despite the quench of the kinetic energy in frustrated systems, localized particles are enabled to move by dissipation, even if this dissipation is strictly local, i.e. each lattice site couples to its own bath. We demonstrate that local Markovian dissipation in a crystal with a flat band permits the transfer of excitations between independent frustrated states and generates coherences between these states. Any mobility due to mixing of flat and dispersive bands can be neglected provided their energy separation is much larger than the systemenvironment couplings. Whilst the lattice sites interact locally with their own baths, frustrated states share these baths such that second-order processes, whereby an excitation hops into the bath and back into the system, allow excitations to tunnel from one frustrated state to another and generate first order coherence between these states. The generation of coherence by dissipation also distinguishes this scenario from pure dephasing which would preserve the number of excitations. Moreover, the dissipation terms we consider here are different from dissipators that can be engineered via coherent driving [33] as they do not conserve the number of excitations. We detail the model in Section II and show that local Markovian baths generate dissipative coupling terms between the flat band states of a frustrated lattice. To explore dissipation-induced mobility we introduce a localised pump and investigate the excitation density in lattice sites away from the pump in the stationary state. In Section III, we find an increased density of excitations that, in the absence of interactions, decays exponentially in the distance from the location of the pump. We also investigate coherence between localised Wannier states and find that first order coherence is generated by the dissipation even for an incoherent input. This coherence is long range but is reduced in the vicinity of incoherent pumps. In Section IV, we show that these results cannot be solely attributed to direct coupling between adjacent sites and that long-range incoherent hopping processes need to be taken into account in order to accurately describe the mobility. We then explore the effect of interactions in Section V, showing that dissipation-induced mobility and coherence is a generic effect for frustrated systems. We state our conclusions in Section VI and discuss possible experimental platforms where these effects can be explored.

II. MODEL
The sawtooth lattice consists of triangles connected to each other at the apex of each triangle, see Fig. 1 where a i (b i ) are the annihilation operators on site i for the A (B) sublattice and with ω 0 being onsite energy of each site (we set = 1 and all lengths are expressed in units of the lattice spacing). The on-site interactions are of the form where U A(B) is the interaction strength for multiple occupation of sites on the A (B) sublattice. The driving Hamiltonian will be specified in Eq. (4) below. We start by considering the non-interacting case U A = U B = 0 and return to investigate the influence of interactions later. The Hamiltonian of Eq. (1) can be written in terms of decoupled Bloch modes with frequencies E k = ω 0 + t cos k ± t 2 cos 2 k + 2t 2 (1 + cos k). In the limit t → √ 2t, the lower-energy band of the Hamiltonian Eq. (1) becomes flat with a gap of 2t to the dispersive band. The kinetic energy in the flat band is quenched due to geometric frustration. This is clearly seen by writing H t in a basis of Wannier states which are exponentially localised around the unit cell i. The tight-binding bosonic operators can be expressed in terms of operators W † i (W i ), that create (annihilate) excitations in the Wannier states, as where w A (r) = √ 2 2π π −π dk cos(k/2)e −ikr e ik/2 / √ cos k + 2 and w B (r) = −1 2π π −π dke −ikr / √ cos k + 2 are independent of t due to the flat-band criterion. If the band gap is sufficiently large, excitations into the dispersive band are negligible and we can project the Hamiltonian onto the flat band [26] such that H t = (ω 0 − 2t) i W † i W i , see App C for details. Here, one clearly sees that the localised Wannier states decouple and there is no propagation in H t .
To explore mobility and coherence in this drivendissipative regime, we first consider a scenario where one Wannier orbital is pumped by a coherent drive with amplitude Ω W and frequency ω D . In a frame rotating at ω D and ignoring rapidly rotating terms, for U A = U B = 0 we get where ∆ = ω 0 − 2t − ω D is the detuning of the drive from the flat band and i = 0 is the pumped site. We consider a drive on resonance with the flat band so ∆ = 0.
In an experiment, this drive of a Wannier state can be implemented by driving the A and B sites coherently with a frequency in resonance with the flat band and with suitable amplitudes such that The dissipation of excitations from the lattice is treated in the Born-Markov approximation with independent baths for each of the lattice sites. The system is thus described by the master equationρ = ρ} is a Lindblad-type dissipator and {·, ·} denotes the anti-commutator. Expressing the dissipation in terms of the Wannier states, we geṫ Whilst the reservoirs of the A and B sites are independent, in the Wannier basis adjacent states share a common reservoir which enables a non-local purely dissipative coupling. In turn, when κ = 1, the non-local dissipation terms γ l for l = 0 vanish and individual Wannier states decouple.

III. NON-INTERACTING REGIME
In the absence of interactions, the master equation Eq. (5) can be solved using Ehrenfest's theorem d Ô /dt = Ôρ = Tr{−iÔ[H, ρ] +ÔD[W ]ρ}. As the solutions of Eq. (5) are Gaussian states, they are entirely specified by the first and second order moments of the Wannier operators W i . Using the non-interacting version of the master equation, the equations of motion for the expectation values for W i and W † i W i+j are given by where δ i,j is the Kronecker delta. Eqs. (7) and (8) provide a set of coupled equations that we solve numerically. The non-local dissipation coefficients γ l decrease exponentially with l so we introduce a cutoff at l > 10 for which we set γ l = 0. Also, we expect long range correlations to be negligible so we make the approximation W † i W i+j = 0 for j > 10. We have checked that the results presented here are converged with respect to these cutoffs. A finite size lattice was used which was sufficiently large for boundary effects to be insignificant.
In order to investigate dissipation-induced mobility in this frustrated lattice, we pump only the site at i = 0. The advantage of pumping in the Wannier basis over the site basis is that it is easier to quantitatively disentangle the effect of the non-local dissipator. Pumping a single site on the sawtooth lattice results in a pump profile of finite width in the Wannier basis which obscures the propagation of the excitations enabled by the dissipation.  i W i normalised by the density on the pumped site on a logarithmic scale. Strikingly, we see that the excitations move away from the driven site even though there are no coherent processes which couple the Wannier states. The transport in this system is enabled by the dissipator, which is achieved solely through coupling of the frustrated sawtooth lattice to local Markovian baths. We will show that this mobility cannot be explained in terms of a diffusion process in Sec IV.
This effect can also be clearly seen in the excitation densities of the original lattice. Fig. 2b shows the density on the B sites b † i b i for the same drive intensity pattern as above. In the lattice site basis, the drive amplitudes Ω A,i and Ω B,i are nonzero for i = 0 as we are pumping a single Wannier state. However, while the excitation density does not vanish on the sites neighbouring the central site when there is no non-local dissipation (κ = 1), it is clear that the density profile changes significantly as the non-local dissipation terms increase. This demonstrates that the mobility of the excitations is indeed due to dissipation-enabled transport and that the effect is not an artifact of the transformation into the Wannier basis.
The normalised density decreases exponentially and can be approximated by N i = exp(−|r i |/ξ) for i = 0. The decay length of the density profile can thus be extracted using adjacent densities ξ i = |log 10 N i − log 10 N i+1 | −1 and is shown in Fig. 3 as a function of κ for i = 4. As κ decreases, the non-local dissipation rates γ l increase which leads to a divergence in ξ as κ → 0. In this limit, the dissipation rates on the B sublattice tend to zero which leads to the formation of the dark state The decay length of the density of excitations in the lattice introduced by the pump diverges as this dark state extends over the entire lattice. In the limit κ → ∞, where the dissipation rate on the A sites tends to zero, the corresponding dark state i (−1) i a † i |0 would involve contributions from the dispersive band and is thus not accessible with our driving mechanism.
B. Long-range first-order coherences The local dissipation in the site basis can also generate first order correlations W † j W l in the sawtooth lattice. Remarkably such coherences even emerge when the coherent input is removed Ω W = 0 and the system is incoherently pumped. To describe the incoherent pumping, we add a term i x=a,b ρ} to the master equation (5) where P x,i is the strength of the incoherent pump. The density profile generated by a local incoherent pump is shown in Fig. 2c and exhibits similar exponentially decaying extension as for a coherent input. Numerical simulations show that the first order coherence g (1) (j, l) = W † j W l / W † j W j W † l W l for coherent driving of the Wannier state W 0 is perfectly correlated with g (1) (j, l) = (−1) |j−l| as the density matrix of the system is a product of coherent states with alternating phases. We thus find long-range first-order coherence despite an exponentially decaying spatial density profile. This indicates that the excitations form a condensate whose extension is the decay length ξ of the density profile.  The same function for a system where the coherent input is replaced with an incoherent drive on the B sublattice at site i = 0 is shown in Fig. 4. The dissipationinduced mobility continues to generate significant coherences even in the absence of a coherent input. As demonstrated for the coherent case, mobile excitations propagate coherently through the lattice. However, an incoherent drive also counteracts the build-up of the first-order coherence g (1) (j, l). As the incoherent pump is strongest for the Wannier site W 0 , the coherence is reduced around g (1) (0, l) as seen in Fig. 4. In contrast to recent experimental observations [18], the coherences in our system are thus longer range, due to the absence of disorder.

IV. MODELLING THE DECAY LENGTH OF THE DISSIPATION-INDUCED MOBILITY
In this section we clarify that the dissipation-induced mobility which we find cannot be attributed to a diffusion process, see Fig. 5a, as this fails to describe the decay length of the density correctly. We continue by examining which of the dissipative transfer processes of Eq. (5) contribute significantly to the mobility by comparing the predicted density distributions for two approximate models to our numerical results.
Considering only the direct transfer of excitations from the pumped site to unpumped remote sites through longrange non-local dissipative incoherent tunnelling events, see Fig. 5b, is not sufficient to explain our findings as it predicts that the decay length ξ is independent of κ. A refined version of this model, whereby a corrected density on the pumped site acts as an effective drive for two remote sites with nearest-neighbour hopping, see Fig. 5c, provides a good approximation for the density decay length ξ for a wide range of κ. The model breaks down as κ → 0 where the contribution of the longer range non-local dissipation processes which we have ignored becomes significant.

A. Single site diffusion model
We begin by attempting to model the density profile by approximating the motion as a diffusive random walk where excitations can hop incoherently between neighbouring sites, see Fig. 5a. The continuity equation for the density on site i for a driven-dissipative system is given by where ∂N i /∂t| L/R is the number of particles on site i moving to the left and right respectively, Γ is the dissipation rate and F (x i ) is a source distribution. The net particle flow to the site i+1 in a time interval ∆τ is given by where J is the rate at which particles hop between adjacent sites and similarly for ∂N i /∂t| L . Therefore, the diffusion equation for Eq. (9) is given by where we note that ∆x 2 /∆τ is some scaled diffusion constant D. The steady-state solution of Eq. (12) without sources is given by N i = Ae −xi/Ξ +Be xi/Ξ with the decay length Ξ = JD/Γ. Assuming that the incoherent dissipative coupling can be modelled as a diffusive process, the most natural values for the dissipation and hopping rates are Γ ∝ γ 0 and J ∝ γ 1 . Therefore, the functional form of the Ξ would be which differs from the ξ obtained from the numerical simulations. If we include the direct non-local dissipative coupling of the pumped site to remote unpumped sites as a source term for Eq. (9), which will have a functional form F (x i ) ∝ e −|xi|/ξΩ where ξ Ω is the decay length of the direct coupling, then this will affect the steady-state density by introducing a particular integral with decay length ξ Ω which will not affect the functional form of Ξ(κ). Indeed, this will prevent the density distribution from having the exponential distribution observed in the simulations.

B. Direct non-local dissipative coupling model
As the dissipation-induced mobility we observe cannot be understood in terms of diffusion through nearestneighbour hopping we will now examine which dissipative transfer processes play an important role. The non-local dissipation is long ranged so a simple approximation is to neglect the dissipative coupling between sites which are not pumped such that transport from the pumped site i = 0 to another site i = j is solely mediated by direct non-local dissipative coupling γ j between the central pumped site and the remote unpumped sites, see Fig. 5b.
We can calculate the density distribution for this model by solving Eqs. (7) and (8) for a hypothetical two-site system consisting of the sites i = 0 and i = j only. The zeroth site is pumped with a drive strength Ω W so these equations give the steady-state density on the jth site where γ 0 γ j due to the exponential localisation of the Wannier states so the decay length is given by In the limit κ → 1, the non-local dissipation rates vanish and the density distribution is predicted well by this model.

C. Effective drive model
The direct coupling model was accurate in the limit κ → 1 so we can hope that a more accurate model will be obtained by incorporating nearest neighbour hopping. Therefore, we solve the equations of motion for the sites i = {0, j, j + 1} where j > 1 assuming that these two sites are not directly coupled to the pumping site with nearest neighbour dissipative coupling terms, see Fig. 5c. We make the assumption that for i = {j, j + 1} the only relevant term in Eq. (7) is from the pumped site l = −j such that l ={0,1} γ l W j+l ≈ γ −j W 0 . Furthermore, we assume that the field on the pumped site is only affected by its nearest neighbours so we can treat W 0 as a constant.
We thus calculate the amplitude of the pumped site W 0 by solving Eq. (7) for the pumped site and its nearest neighbours and use this result to calculate the densities on sites i = j and i = j + 1. Details can be found in App. D and result in a decay length of As we see in Fig. 3, this approximation works well for a wide range of non-local dissipation strengths. However, as κ → 0, the approximation breaks down and the decay length is underestimated due to multiple hopping processes which have not been included. We thus conclude that for κ ≈ 1, the direct process from the pumped site to the considered site plays the most significant role. As κ decreases, the transfer of particles from the neighbouring site to the considered site becomes relevant. The smaller κ becomes, the more these indirect processes play a role and the more sites must be included to provide an accurate approximation of the density distribution.

V. STRONGLY-INTERACTING REGIME
After investigating the non-interacting case, we now turn to explore the effect of interactions on dissipationinduced mobility and coherence. As with the dissipator, the projection of the interaction Hamiltonian H U onto the Wannier basis results in non-local interaction terms, see App. E 1 for a full derivation. These non-local interactions can also cause propagation, as has been explored in detail [26,[29][30][31][32].
The interaction induced terms which cause coherent propagation are least effective in the limit of very strong interaction, where propagation from terms such as (W † i ) 2 (W i+1 ) 2 is suppressed. To explore the effects of the dissipation, we thus consider this limit, where the dominant on-site interactions suppress multiple Wannier state excitations, and for simplicity set U A = 0. Truncating the system to the single-excitation subspace, the Hamiltonian in Eq. (2) in the Wannier basis can be approximated by The leading-order contributions of H U are the onsite interaction, nearest-neighbour cross-Kerr interactions, which cannot move an excitation from one site to another, and density-assisted tunnelling, whereby an excitation can coherently tunnel if there is an excitation on a nearby site. The density-assisted tunnelling terms contribute to the mobility for our driven-dissipative sawtooth lattice and could obscure or even dominate over the contribution induced by the non-local dissipator. We now show that the joint effect of the interaction terms suppresses the mobility in the strong interaction regime and that, even here, the main contribution to the mobility originates from the dissipation. To this end, we calculated the density distribution in the strongly interacting regime using a variational matrix product operator (MPO) method [34][35][36] (see App. E 4). The interactions modify the density distribution which, in contrast to the non-interacting case, no longer decays exponentially away from the pumped site, see Fig. 6a. To study how H U affects the dissipation-induced mobility, we thus use the fraction of excitations in the unpumped Wannier states f = ( i =0 N i )/( i N i ) as the relevant figure of merit (f = 0 for uncoupled Wannier states). Fig. 6b shows this non-local density fraction as a function of κ for a range of U 1 for which the truncation to the single-excitation subspace is valid (see App. E 2). For fixed κ, the mobility of the excitations decreases with increasing interaction strength. We attribute this reduction in the mobility to the cross-Kerr interactions, which shift the energy of the neighbouring site, detuning the non-local dissipative transition and preventing excitations from moving between sites. The effect of the density-assisted tunnelling terms is insignificant as setting U 2 = 0 leads to a variation of f of the order of 10 −2 . Moreover, the dependence of f on U 1 is much weaker than its dependence on κ, showing that the non-local dissipator is the dominant contributor to the steady-state of the  system in the limit of strong interactions. We also computed the coherences between lattice sites for U 1 = 100, Ω W = γ A and κ = 0.1, see Fig. 7. Similarly to the incoherent drive, interactions counteract the build-up of first order coherence. Since the density of excitations is highest for the pumped site, interactions are most effective on this site, which reduces g (1) (j, l) around j = 0. Our results demonstrate that the generation of coherences due to dissipation-induced mobility is a general phenomenon and not restricted to non-interacting systems.

VI. CONCLUSIONS
We have shown that local Markovian dissipation can induce mobility and long range coherence in frustrated lattice systems in the absence of kinetic energy. An experimental demonstration of the effect we predict here could be realised in any driven-dissipative lattice with a flat band. For example, in photonic waveguides lattices [15], the couplings between sites can be engineered to realise a flat band system where defects and disorder can be kept low enough to not affect the dynamics of the lattice excitations. Neighbouring waveguides can be used to simulate the local Markovian baths with a strong degree of control over the environmental dynamics as has recently been demonstrated [37]. Additionally, microwave resonator arrays in superconducting platforms with their exquisite precision in fabrication and control provide a suitable system. An imbalance in the dissipation rates of A and B sublattices is expected, or can be engineered, in many situations and could be made stronger via Purcell enhancements. Whilst we have discussed steady-state density distributions, the effects we find could also be studied dynamically, where varying the imbalance in dissipation rates would change the rate at which a localised excitation disperses across the lattice. Photonic waveguide lattices [15] would be particularly well suited for such experiments. The equations of motion for the driven-dissipative sawtooth lattice is given by the master equationρ = ρ} describes Markovian loss of an excitation from site i. For simplicity, we will assume that the losses from all the A sites and all B sites are the same but the loss rates for A and B sites are not necessarily equal ie. γ A,i = γ A and γ B,i = γ B .
If we take Eq. (3) and plug it into the dissipator where we note that w A (r) is real as it is a Fourier transform of a symmetric function. Similarly, for the dissipa- Combining these two equations, we have where γ eff j,l = γ eff l only depends on the separation l of the Wannier states and cos k(l + 1) + cos k(l − 1) cos k + 2 dk If γ A and γ B are equal, we can take the sum over i which gives a local dissipator due to the Wannier state orthogonality relation (A11) and all non-local dissipation terms disappear. When the dissipation rates are not equal, we can evaluate the integral f l using the residue theorem. For simplicity, let us assume that l ≥ 0 and we note that f l = f −l . First, we make the standard substitution z = e ik such that π −π cos kl cos k + 2 dk = C z 2l + 1 where the integration contour C runs anticlockwise around the unit circle. This integral has three poles: two simple poles at z = −2 ± √ 3 and a pole of order l at z = 0. The pole at z = − √ 3 − 2 lies outside the unit circle so it does not contribute to the integral. The residue at z = √ 3 − 2 is calculated using the standard method for a simple pole For the pole at z = 0, we use the Laurent expansion of the integrand and take the prefactor of the z −1 term. The lowest order term in the Laurent expansion when z 2l is in the numerator is of order z l so we only need Therefore, the integral in Eq. (A12) can be evaluated using the residue theorem to give π −π cos kl cos k + 2 dk = 2π where the absolute value of l is taken as we know that the integral is symmetric under the transformation l → −l.
For reference, some numerically integrated values of f l are given below: We see that the dissipator contains non-local terms but that the dissipation rate for coupling to adjacent sites drops off rapidly. There are negative dissipation rates but the dissipation rate matrix γ eff j,l is positive. To connect to the non-local dissipation coefficients for the Lieb lattice calculated in the following section of this appendix, we also note that f l can be expressed in terms of the regularised hypergeometric function f l = 3F2 (1/2, 1, 1; 1 + l, 1 − l; −2)

Appendix B: Wannier states of the Lieb Lattice
In the main text, we have focused on the sawtooth lattice, which is the simplest one-dimensional frustrated system. However, the conclusions we derive in this paper are equally applicable to other frustrated lattices where the master equation for the Wannier states of the flat band have different dissipation rates for the different sublattices, as in Eq. (5). Notably, the Lieb lattice, studied in some contemporary works on driven-dissipative frustrated systems [18,24,25], exhibits the same behaviour.
The Lieb lattice consists of a one-dimensional chain of harmonic oscillators with nearest-neighbour coupling J decorated with a coupling g to an additional harmonic oscillator for every other site in the chain. The Hamiltonian for this system is given by where the unit cell is labelled with A, B and C sites, ω A , ω B and ω C are their respective onsite energies and a i , b i and c i are the annihilation operators for the respective sublattices of the unit cell i. Transforming Eq. (B1) into the momentum basis and expressing distances in units of the lattice spacing one gets For the flat band, the normalised eigenstates can be calculated using the ansatz ψ † k = cos θ k b † k + e iφ k sin θ k c † k from which we find that φ k = k/2 and tan θ k = −(g/2J cos(k/2)) −1 . The Wannier coefficients are given by e ikr e ik/2 1 + g 2 /4J 2 cos 2 (k/2) (B4) Therefore, as with the sawtooth lattice, the creation operators of the Lieb lattice can be expressed as These are eigenstates of the Lieb Hamiltonian Eq. (B1) so projecting onto the flat band we have H = j ω B W † j W j . If we introduce local, independent dissipation on the B and C sublattice, we obtain a similar non-local dissipator to the one used for the sawtooth lattice. The dissipation rate of the A sublattice is irrelevant as the flat band eigenstates have no weight of the A sites.
The derivation of the non-local dissipation coefficients follows from Eqs. (B5) and (B6) in exactly the same way as for the sawtooth lattice. Therefore, where γ B and γ C are the dissipation rates on the B and C sites respectively and κ = γ B /γ C . We can evaluate the terms in this sum i e i(k+k )ri e −ik rj 1 + (2J cos(k/2)/g) 2 1 + (2J cos(k /2)/g) 2 dkdk (B8) = 1 2π π −π cos(kj) 1 + (2J cos(k/2)/g) 2 dk (B9) i e i(k+k )ri e −ik rj e i(k+k )/2 1 + (g/2J cos(k/2)) 2 1 + (g/2J cos(k /2)) 2 dkdk (B11) where a = (2J/g) 2 . We can express the non-local dissipation coefficients in a similar form as for the sawtooth lattice by noting that where (x) k is the Pochhammer symbol and we have used the identity (x) k = Γ(x + k)/Γ(x), which implies that (x + 1) k−1 = (x) k /x and (1) k−1 = (1) k /k. Additionally, we obtain the delta-function 1/Γ(1+j)Γ(1−j) → δ j,0 as j takes only integer values, Γ(1) = 1 and |Γ(1 − |j|)| → ∞ for all j = 0. Therefore where the geometric factor f j (a) = 3F2 (1/2, 1, 1; 1+j, 1− j; −a) is now a function of the tunnelling rates J and g. Fig. 8 shows the distance dependence of the magnitude of the non-local dissipation rate for a range of values of a = (2J/g) 2 . We see that as J/g increases, the characteristic length scale of the incoherent transport increases. This can be understood intuitively as this coupling increases the kinetic energy along the undecorated lattice (i.e. for the system where the coupling to the C sublattice vanishes) so excitations can more easily flow along the lattice, even through dissipative processes.

Appendix C: Derivation of the Wannier state master equation
In this appendix, we give a derivation of the Lindblad master equation Eq. (5) in the Wannier basis from an initial, general system-environment Hamiltonian and show that it is valid in the limit of weak system-environment coupling.
Initially, let us assume that there is no driving potential Ω = ω D = 0. We assume that the system consists of the original sawtooth lattice where each site is independently coupled to a large set of harmonic oscillators which act as independent baths for each of the lattice sites. The total Hamiltonian for this setup is given by with where g a k,i is the coupling strength of the ith A site to the harmonic mode k, c a k,i is the creation operator for this mode and ω a k,i is the mode energy and respectively for the B sites.
To find the master equation for the Wannier states, we transform to the Wannier basis using Eq. (3) and trace out the environmental degrees of freedom neglecting the dispersive band. For the Wannier states, H S and H I are defined as where, for clarity, we have abbreviated the Wannier state coefficients w A (r i − r j ) = w ij A and w B (r i − r j ) = w ij B . To calculate the master equation in the Born-Markov approximation, we start from [38] where H I (t) = e i(H S +H E )t H I e −i(H S +H E )t is the interaction Hamiltonian in the interaction picture, ρ(t) and ρ E is the reduced density matrix of the system and environment respectively and Tr E {·} is the trace over the environmental degrees of freedom. As both H S and H E are diagonal, it is simple to show that The dispersive band couples similarly to the bath through a separate interaction Hamiltonian where ∆ q , V q and v iq A(B) are the energies, annihilation operators and matrix coefficients for the Bloch states of the dispersive band respectively. By virtue of the orthogonality relations of the bath operators in Eq. C13-C15, the terms in Eq. C7 which couple the flat and dispersive bands by including both H I and H V I only contain contributions where the flat and the dispersive band couple to the same bath mode. Since the energy difference between the flat and dispersive band largely exceeds the strength of the system bath couplings, these terms can be neglected in a rotating wave approximation. As a consequence, it is a good approximation to restrict the description to the flat band interacting with the environments as we do in our approach. This justifies our assumption that the dispersive band does not affect the mobility of excitations in the flat band.
Evaluating the trace for each term of the double commutator is tedious so here we just focus on the first term where the trace over the environmental harmonic modes, assuming that they are in thermal states, allows us to use the orthogonality relations Here n a k,i is the occupation of the kth mode of the bath of the harmonic oscillator coupled to the A site of lattice site i and similarly for the B sites and we have defined The resulting double commutator is given by Separating F 1 jj and F 2 jj into real and imaginary parts allows us to complete the master equation where D[A, B]ρ ≡ 2AρB − {BA, ρ} and The master equation Eq. (C19) also contains a Lamb shift with We will now make a couple of approximations about the environment and the system's coupling to it. Firstly, we assume that the environment is in its ground state and contains no excitations, i.e. n a k,i ≈ 0 and n b k,i ≈ 0 so F 2 j,j ≈ 0. Whereas generalizations to finite temperature environments would also yield the mobility we find, the assumption of a zero temperature environment is a very good approximation for systems hosting optical or microwave photons. Secondly, let us assume that the spectrum and the coupling strength of each of the sites to its local environment is the same. This means that g a k,i , g b k,i , ω a k,i and ω b k,i are all independent of i. However, we allow these properties to be different between the A and B sites such that g a k,i = g b k,i and ω a k,i = ω b k,i so that the dissipation rates on the A and B sites can differ. The non-local dissipation rates can be calculated directly from this calculation 3 is calculated using the residue theorem as shown in App. A.
An inhomogeneity in the Lamb shifts described in Eq. (C22) would lead to a kinetic energy term an thus to additional mobility. These Lamb shifts are equivalent to local frequency shifts, and are therefore irrelevant for most experimental situations.
In photonic waveguide lattices, the propagation distance along the waveguide takes on the role of time and the propagation constant the role of the resonance frequency. The latter is set by the real part of the refractive index which is unaffected by the loss channels that are patterned into the chip via additional waveguides. For lattices of coupled resonators, in turn, one aims to build mutually resonant lattice sites. Therefore the lattice sites with larger damping rates should be built to have a correspondingly adjusted frequency, which is not more difficult than building them at exactly equal frequencies. Any residual frequency disorder can in many experimental situations be compensated by tuning the individual lattice sites, e.g. via ac-Stark shifts. Since any measurement of the resonance frequency of a resonator always includes the Lamb shift, disorder in these shifts is then automatically compensated for when tuning the device.

Appendix D: Decay Rate Derivation for the Effective Drive Model
In this appendix, we derive the decay length of the density profile for the non-interacting system using an effective drive model. This model involves solving the equations of motion for the three site system i = {0, j, j + 1} where j > 1 assuming that the two unpumped sites do not affect the field amplitude of site i = 0, which is a constant determined by the three site system i = {−1, 0, 1}. This model combines nearest-neighbour and direct non-local dissipative coupling terms and results in an accurate prediction of the decay length ξ for a large range of dissipation rate asymmetry κ.
To calculate the field amplitude for the pumped site, we start by solving Eq. (7) for the system i = {−1, 0, 1}. The equations of motion are given by where we have exploited the parity symmetry around i = 0 and the relation γ l = γ −l . Solving Eqs. (D1) and (D2) for the pumped site's steady state amplitude gives which we can insert into the steady state of Eq. (7) along with the assumption that l ={0,1} γ l W j+l ≈ γ −j W 0 which we made in Sec.
We see that the pumped site acts as an effective drive for the two site system with These effective drives can be used in the equations of motion for the correlations of the Wannier state operators in the steady state, namely We can solve Eqs. (D4)-(D11) to find approximate densities on the sites j and j + 1 from which we can calculate the decay length As discussed in Sec. IV C, this provides a good approximation for the exact density decay length ξ for a large range of κ. U A n a † n a † n a n a n =U A i,j,l,m,n dkdk dqdq · cos(k/2) cos(k /2) cos(q/2) cos(q /2) n e −irn(k+k +q+q ) (cos k + 2)(cos k + 2)(cos q + 2)(cos q + 2) · e ikri e ik rj e iqr l e iq rm e i(k+k +q+q ) dk dqdq cos((k + q + q )/2) cos(k /2) cos(q/2) cos(q /2) (cos(k + q + q ) + 2)(cos k + 2)(cos q + 2)(cos q + 2) · e ik (rj −ri) e iq(r l −ri) e iq (rm−ri) Π (k + q + q ) 2π and for the B sites e ikri e ik rj e iqr l e iq rm n e −irn(k+k +q+q ) (cos k + 2)(cos k + 2)(cos q + 2)(cos q + 2) dk dqdq e ik (rj −ri) e iq(r l −ri) e iq (rm−ri) (cos(k + q + q ) + 2)(cos k + 2)(cos q + 2)(cos q + 2) dk dqdq cos((k + q + q )/2) cos(k /2) cos(q/2) cos(q /2) (cos(k + q + q ) + 2)(cos k + 2)(cos q + 2)(cos q + 2) · e ik j e iql e iq m · Π (k + q + q ) 2π (E10) dk dqdq e ik j e iql e iq m (cos(k + q + q ) + 2)(cos k + 2)(cos q + 2)(cos q + 2) are the effective interaction strengths and Π(x) is the tophat function of unit width centered at zero. The U eff,A j ,l ,m and U eff,B j ,l ,m are symmetric on interchange of the indices j , l and m and U eff,B j ,l ,m is also symmetric for j → −j , l → −l and m → −m . These effective interaction strengths decrease exponentially with increasing j , l and m as they consist of sums over the Wannier coefficients, which are exponentially localised.

Truncating to the single excitation subspace for MPO simulations
The sum in Eqs. (E5) and (E9) describe an infinite number of interaction processes. However, these interactions can be divided into three types: density-assisted tunnelling, cotunnelling and density-density interactions. Density-assisted tunnelling terms produce coherent coupling between two sites which is activated when there is an excitation either on one of the tunnelling sites or on a neighbouring site. An example is given by j = 1, l = 1, m = −1 which describes coherent tunnelling between the sites i and i − 1 which occurs only when there is an excitation on site i + 1. Cotunnelling terms describe two excitations on different lattice sites tunnelling together onto two other lattice sites. An example would be when j = −1, l = 1, m = 2 which corresponds to a tunnelling of an excitation on the i and i − 1 sites to the i + 1 and i + 2 sites. Finally, densitydensity interactions terms with a cross-Kerr non-linearity such as for j = 1, l = 0, m = 1 which gives an energy penalty for having excitations on neighbouring sites.
Whilst the density-density interactions cannot transfer excitations between sites, the density-assisted tunnelling and cotunnelling terms both introduce coherent transfer processes. We want to investigate how interactions affect the dissipation-induced mobility we examine in this paper. We expect that the cross-Kerr non-linearity will decrease the mobility as it introduces an energy shift on sites neighbouring an excitation. This will detune the neighbouring site and prevent excitations moving into these states due to energy conservation. On the other hand, the density-assisted tunnelling and cotunnelling terms should increase the mobility.
This proliferation of possible transfer processes makes it difficult to determine which processes are affecting the density distribution when including all interaction terms in Eqs. (E5) and (E9). In order to reduce the complexity of this problem, we choose to work in the stronglyinteracting limit where the probability that a site contains two or more excitations is vanishingly small. This means that the operator H U vanishes when two Wannier state annihiliation or creation operators operate on the same site.
This truncation also facilitates our calculations of the density profiles in the interacting regime using a variational Matrix Product Operator (MPO) method [34][35][36]. The computational resources necessary for this method decrease significantly if the system can be described by a reduced Hilbert space on each lattice site. For our truncation approximation to be valid, we need U 0 to be large enough to prevent the pump doubly-exciting the site at i = 0.
To see when this approximation is valid, we calculate the second order correlation function for the pumped site W † 0 W † 0 W 0 W 0 assuming that all coherent and incoherent processes coupling it to other sites vanish. This gives us an upper bound on the probability that the pumped site contains at least two excitations as the addition of non-local terms will only decrease the density on the pumped site. For this case, W † 0 W † 0 W 0 W 0 can be calculated exactly using an approach by Drummond and Walls [39,40]. The result is plotted in Fig. 9 as a function of U 0 /Ω W,0 and γ 0 /Ω W,0 .
Using the simulations performed in the main text allows us to perform a consistency check on the truncated Hilbert-space approximation for the parameters we have used. The excitation density of the site next the pump for the non-interacting system is observed to be in the range of 10 −2 − 10 −3 for κ > 0.1. We want the probability of finding a double-excitation to be smaller than this. In our simulations, we set γ A = 1 and Ω W,0 = 1 so for W † 0 W † 0 W 0 W 0 10 −4 , Fig. 9 shows that we require U 0 2.0. This provides us with a minimum onsite interaction strength U 0 which we must use to guarantee the validity of the single excitation subspace approximation. Importantly, the non-local interaction strengths where j , l and m are not equal to zero are related to U 0 so we are unable to make them arbitrarily small. However, we can calculate minimum value for these non-local terms which correspond to a regime where the truncation to the single excitation subspace is valid.

Approximate Interaction Hamiltonian for the Single Excitation Subspace
As described above, in our MPO simulations we need to truncate the Hilbert space of each lattice site to the subspace of at most one excitation. Therefore, we exclude terms with two annihilation or creation operators on the same site, with the exception of the onsite interaction. where {j , l , m } consists of a set of terms grouped using the symmetries of U eff,A j ,l ,m and U eff,B j ,l ,m described in the previous section of this appendix. Terms which bring the state out of the single excitation subspace are not included.
Including interaction effects for both A and B sublattices is an unnecessary complication, as it only leads to minor modifications in the weighting of the interaction strengths. We have tested that these small modifications only have a negligible effect. We choose U A = 0 and U B = 0 as this leads to the strongest possible coherent propagation terms and allows us to show that even here the incoherent mobility is the dominant effect.
In the strongly-interacting regime, the leading order terms in Eq. (E1) in the Wannier basis is the onsite energy with {j , l , m } = {0, 0, 0}, the nearest-neighbour cross-Kerr interaction with {j , l , m } = {1, 1, 0}, the next-to-nearest-neighbour cross-Kerr interaction with {j , l , m } = {2, 2, 0} and the density-assisted tunnelling term with {j , l , m } = {−1, 0, 1}. However, the symmetries of U eff,B j ,l ,m provide multiplicative factors for equivalent terms in the sum in Eq. (E9). The onsite term does not have a multiplicative factor. The cross-Kerr interactions have a fourfold multiplicity, as can be seen for density-density interactions between sites i = 0 and i = 1 (i.e. W † 0 W 0 W † 1 W 1 ) where the following terms all contribute to the interaction strength: For the density-assisted tunnelling terms in the symmetry set {j , l , m } = {−1, 0, 1}, where the density operator is on the central site and the excitation hops over it, e.g. W † −1 (W † 0 W 0 )W 1 where i = 0, j = −1, l = 0 and m = 1. This process has a twofold multiplicity in the sum in Eq. (E9) due to the symmetry j ↔ m .
Considering these leading order terms, we can approximate the interaction Hamiltonian by where U 0 ≈ 0.192U B , U 1 ≈ 0.133U B , U 2 ≈ 0.054U B and U 3 ≈ 0.023U B , which is the approximation to the interaction Hamiltonian H U which we use for the MPO calculations. Using the minimum value for U 0 needed to apply the single excitation subspace truncation which we calculated in the previous section of this appendix, we see that the approximation in Eq. (E13) is valid for U B 10.0.

Agreement of results for interacting and noninteracting cases
The MPO method expresses the vectorised density matrixρ as a matrix product state (MPS) and the Lindbladian L as an MPO. The MPS for this system had a virtual dimension of 40 and the non-local dissipation was cut off at l > 3. To find the stationary state we seek the ground state of L † L in the manner of [36].
We note that for low enough densities, our numerical MPO calculations are in agreement with the results for the noninteracting case. However, a direct comparison of interacting and noninteracting scenarios for higher densities would require the simulation of more states per lattice site which is beyond the capabilities of our MPO simulations.