Destruction of Localization by Thermal Inclusions: Anomalous Transport and Griﬀiths Effects in the Anderson and André-Aubry-Harper Models

We discuss and compare two recently proposed toy models for anomalous transport and Griﬃths eﬀects in random systems near the Many-Body Localization transitions: the random dephasing model, which adds thermal inclusions in an Anderson Insulator as local Markovian dephasing channels that heat up the system, and the random Gaussian Orthogonal Ensemble (GOE) approach which models them in terms of ensembles of random regular graphs. For these two settings we discuss and compare transport and dissipative properties and their statistics. We show that both types of dissipation lead to similar Griﬃths-like phenomenology, with the GOE bath being less eﬀective in thermalising the system due to its ﬁnite bandwidth. We then extend these models to the case of a quasi-periodic potential as described by the André-Aubry-Harper model coupled to random thermal inclusions, that we show to display, for large strength of the quasiperiodic potential, a similar phenomenology to the one of the purely random case. In particular, we show the emergence of subdiﬀusive transport and broad statistics of the local density of states, suggestive of Griﬃths like eﬀects arising from the interplay between quasiperiodic localization and random coupling to the baths


Introduction
Transport properties of quantum many body systems are known to be extremely sensitive to inhomogeneities and spatial disorder. A well understood example is provided by the phenomenon of Anderson localization [1], where diffusion is suppressed by quantum interference above a critical disorder strength, which turns out to be vanishingly small in low dimensions [2]. Another notable example is provided by the case of a quasi-periodic potential incommensurate with the lattice, as in the model of André-Aubry-Harper (AAH) [3,4] where in one dimension a transition from delocalized ballistic transport to localized behavior arises as a function of the potential strength.
An appealing phenomenological interpretation of these anomalies has been proposed, for disordered models, in terms of the existence around the MBL transition of a quantum Griffiths phase [36]. This is characterized by rare inclusions of the insulating (conducting) phase with an anomalously small (large) localization length. Several phenomenological proposals have been put forward to describe this physics, from purely classical resistor-capacitor models with power-law distributed resistances [23,37] to stochastic models for merging thermal and insulating regions which are motivated by strong-disorder renormalization group arguments [38][39][40][41]. Crucially, this Griffiths physics is at the core of our understanding of both the critical properties of the MBL transition and the mechanism for its destabilization through quantum avalanches [42][43][44][45][46][47][48].
In the quasi-periodic case the understanding of transport and thermalization across the MBL transition is even more preliminary. Anomalous transport has been reported in the past years for the interacting AAH model, including claims of subdiffusion and superdiffusion [49][50][51], which have been, however, recently questioned [52]. The existence of a subdiffusive region, prior to the fully localized phase, is an intriguing possibility and its relation with a possible Griffiths phase is a priori not obvious. In fact, the quasi-periodicity is a structured pattern, whereas disordered potentials are random and uncorrelated processes. A better understanding of the robustness of quasi-periodic localized phases against thermal inclusions and their transport properties is therefore urgent.
Unfortunately, the study of large disordered or quasi-periodic interacting many body systems with numerically exact methods remains extremely challenging. It is therefore desirable to find models which capture some of the key physical facets of many-body localizing systems, while being computationally tractable even for large system sizes.
Recently, two toy models have been proposed to describe, at a microscopic level, the interplay between localization and the so-called thermal inclusions, which in truly many-body systems are produced by internal interactions, in the context of MBL with random disorder. The first one, called random dephasing model [53], describes the environment in terms of Markovian dephasing processes able to heat up the system and which are randomly coupled to each lattice site with a given probability. The second one describes thermal inclusions as the coupling to a GOE random-matrix bath implemented through an ensemble of Random-Regular-Graphs (RRG) [54]. In addition to their role for the understanding of MBL physics these models also present an intrinsic interest, as simple settings where questions related to stability of localization with respect to dissipation can be answered in some detail.
It is important to mention here that some recent works have even questioned the existence of a truly MBL phase in the thermodynamic limit [55][56][57][58][59]. Furthermore, several recent numerical results indicate that, even if a genuine MBL transition exists, it is affected by strong finite-size effects, with the transition point moving to larger and larger values of the disorder as larger system sizes are considered [60,61], and eventually following far outside the numerically-accessible crossover between the finite-size MBL regime and thermalization. Yet, the subdiffusive crossover region preceding the transition seem to be very robust and much less affected by finite size effects. In this respect Refs. [62][63][64][65] put forward the idea that that the MBL transition and the anomalously slow sub-diffusive crossover phase preceding it may be driven by distinctly different physical mechanisms. In this paper we focus on this latter regime. In this regard, the advantage of working with simplified toy models is twofold. On the one hand, they allow one to inspect the physics on much larger scales (on which Griffiths effects are believed to be important) compared to the interacting models. On the other hand, they allow one to study the sub-diffusive regime disentangling it from the MBL transition.
In the first part of this work we revisit and compare in detail these two models of thermal inclusions in the context of a simple Anderson insulator. In particular, we compute transport and spectral properties, going beyond the results presented in [53,54]. We then extend these toy models to the quasi-periodic case, with the goal of shedding light on the mechanism for the destruction of their localized phases by thermal inclusions and assessing its degree of universality.
The paper is structured as follows. In Sec. 2 we give an overview of the main results of our paper. In Sec. 3 we define the Anderson and AAH models, the dissipative settings we consider, the main quantities of interest and we also briefly mention how to compute them. In Sec. 4 we present our results for the random disordered case, while in Sec. 5 we discuss the quasi-periodic one. In Sec. 6 we discuss and compare the behavior of the two models and we discuss them in the light of possible connections with the MBL problem. Finally, in Sec. 7 we present our conclusions. We detail in the Appendix the derivations and computations of our analytic results.

Summary of Main Results
We start with an overview of the main results obtained in this work, which will be discussed in detail in Sections 4 and 5.
We consider two models for single particle localization, namely one dimensional quantum particles in presence of (uncorrelated) random or quasiperiodic disorder leading to the Anderson and AAH models defined in Sec. 3.1, and we study the effect of thermal inclusions on their transport and spectral properties. As we discuss in Sec. 3.2 for each of these two models we consider and compare two types dissipative environments: (i) a Markovian Dephasing Bath (MDB) described by the Lindblad master equation with randomly distributed dephasing jump operators, similar to what was done in Ref. [53], and (ii) a GOE bath described by a local coupling to independent Random-Regular-Graphs, as done in Ref. [54]. The advantage of these two settings is that in both cases the effect of the coupling to the bath can be treated exactly using Lindblad equations of motions or the cavity method, respectively, and these lead to numerically exact results for single-particle observables such as the particle current or the local density of states, as discussed in Sec. 3.3.
We first discuss in Sec. 4 the case of the Anderson model coupled to MDB and GOE baths, for which we revisit and complement the results of Refs. [53,54]. In particular, we discuss transport properties such as the scaling of the typical resistivity with system size, a study not performed in Ref. [54], and the statistics of the Local Density of States (LDoS), missing in Ref. [53], thus obtaining a complete picture of the effect of thermal inclusions on the Anderson localized phase leading to a robust sub-diffusive phase. Our results are summarized in the dynamical phase diagram shown in Fig. 5. This study lets us demonstrate that the two ways of introducing dissipation, MDB and GOE, lead to the same physical behavior. Already from this analysis we can conclude that the GOE bath is less effective in thermalizing the system and inducing diffusion.
In Sec. 5 we extend our analysis to the quasiperiodic case, i.e. the AAH model. We discuss how dissipation in the two settings affects its transport and spectral properties. Our results show that in both cases the localized phase of the AAH model turns into a sub-diffusive regime which appears to be much broader for the GOE dissipation than for the MDB one, as shown in the dynamical phase diagram of Fig. 9. This subdiffusive scaling appears also in the statistics of the LDoS.
A comment is in order here. For a truly many-body model in a quasiperiodic potential the spatial structure of the thermal inclusions is not random and is arguably correlated with the potential itself. However, predicting the position of these thermal spots is an extremely hard task, which is essentially equivalent to solving the many-body problem. For this reason in the toy models we work with the dephasing dynamics is assumed to apply on a randomly distributed set of sites, independently of the local potential, both in the Anderson and in the AAH case. One could then argue that, since in our simplified setting the insulating regions and the thermal inclusions are essentially put by hand at random, one cannot draw any conclusion on the universality of the Griffiths picture. Yet, on the one hand one could speculate that manybody interactions, together with the randomness of the initial configuration, give rise to some sort of configurational disorder, as one could see for example by treating those interactions at the Hartree-Fock (HF) level [66,67]. Besides, we would like to point out that the fact that the Anderson model with uncorrelated disorder and the AAH model in the localized regime respond in a very similar way when coupled to a thermalizing system is still interesting and informative. In fact, as explained below, in the Anderson model the formation of rare resonances is a quite "simple" local process: resonances are formed on the sites on which the local disorder is small and the dephasing dynamics is active. Instead the formation of resonances in the AAH case cannot be predicted easily and is likely to involve more complex and non-local processes. Still, the phenomenology of the two models when coupled to the thermal inclusions is essentially the same, for the two kinds of dissipative enviroments that we consider. We believe that this gives some strong indication on the fact that the subdiffusive regime is a very robust feature which appears whenever Anderson localized edgestates are perturbed by thermal inclusions, independently of the details of the microscopic modelization of the bath.

Models and Methods
In this Section we define the models we work with. We first introduce the two non-interacting one-dimensional chains, namely the Anderson tight-binding model with quenched randomness and the André-Aubry-Harper model with a quasi-periodic potential. Next, we explain the two ways in which we introduce the thermal inclusions due to the coupling to an effective environment. Finally, we identify the observables we employ and we briefly describe how we implement them (a thorough analysis is detailed in the Appendices).

Anderson and André-Aubry-Harper Models
Let us consider a one-dimensional tight-binding model of non-interacting spin-less fermions moving on a chain of L sites. The system's Hamiltonian iŝ In Eq. (1)d i (d † i ) are the annihilation (creation) operators acting on the i-th site, the oneparticle Hamiltonian is an L × L Hermitian and tridiagonal matrix, and open boundary condition are applied. The constant hopping t sets the unit of energy, and the inhomogeneous potential h i is taken to be either random and uniformly distributed or quasi-periodic (i.e. incommensurate with the lattice) with = ( 5 − 1)/2 (known as the golden ratio) 1 and λ ≥ 0. These are the Anderson [1] and the André-Aubry-Harper quasi-periodic [4] models, respectively.

Random Coupling to Baths
With the aim of mimicking the effect of thermal inclusions, we consider a random coupling to baths, modeled in two ways which we specify in the following two subsections. The couplings allow for energy relaxation, with a rate γ i that is a random variable distributed according to where p ∈ [0, 1] is a control parameter, and γ sets the dephasing strength (see the cartoons in Fig. 1 and Fig. 2). In particular, the limits p → 1 corresponds to no dephasing at all, -p → 0 is a constant and uniform dephasing acting on every site.
Given a spin chain of length L, (1−p)L is the average number of sites under dephasing. Finally, a reason for setting these coupling terms randomly across the set-up is to mimic the spatial distribution of thermal inclusions in many-body localized systems (at least when considering models with a random potential). For deterministic models the position of a thermal inclusion is not random, although it remains unpredictable. For this reason we keep the random coupling distribution (4) even when focusing on André-Aubry-Harper chains.

Markovian Dephasing Bath
In order to include the effect of the environment while keeping the system tractable we work in the framework of open Markovian quantum systems described by the Lindblad equation for the system density matrix [68], i.e.
where the first term represents the coherent evolution of the density matrix due to the system's HamiltonianĤ, while the second term accounts for incoherent processes due to the coupling to the environment and it is described by a dissipator D[ρ t ] acting on the density matrix. In this work we consider three types of dissipative processes, which give rise to a dissipator of the form egions of strong disorder result in bottlenecks ransport. Also quasiperiodic system exhibits transport, but the correlated nature of the ntial renders the Griffiths phase explanation cated guess. The stability of a subdiffusive r to the fully localized phase, is a fascinating supported by recent numerical works [22,23]. able to find a toy model which capture the key pect of many-body localizing systems, while able. In these notes, we consider the ideas in the non-interacting Aubry-André-Harper el. Specifically, the system is evolved coupled undary driving, and to percolated on-site deach site is subject to dephasing, except for issa.it with probability p, 1 p. The limiting case p = 0 correspond to dephasing at every site, while p = 1 to no dephasing.[To Do: larger fonts]

II. MODELS
The general evolution of open quantum systems coupled to a Markovian environment is encoded in the Lindblad equation order result in bottlenecks siperiodic system exhibits correlated nature of the riffiths phase explanation stability of a subdiffusive zed phase, is a fascinating numerical works [22,23]. del which capture the key localizing systems, while es, we consider the ideas ting Aubry-André-Harper system is evolved coupled to percolated on-site deto dephasing, except for with probability p, 1 p. The limiting case p = 0 correspond to dephasing at every site, while p = 1 to no dephasing.[To Do: larger fonts]

II. MODELS
The general evolution of open quantum systems coupled to a Markovian environment is encoded in the Lindblad equation te potential renders the Griffiths phase explanation n educated guess. The stability of a subdiffusive n, prior to the fully localized phase, is a fascinating bility, supported by recent numerical works [22,23]. is desirable to find a toy model which capture the key ical aspect of many-body localizing systems, while treatable. In these notes, we consider the ideas f. [24] in the non-interacting Aubry-André-Harper ) model. Specifically, the system is evolved coupled to a boundary driving, and to percolated on-site deing. (Each site is subject to dephasing, except for rkesh@sissa.it

II. MODELS
The general evolution of open quantum systems coupled to a Markovian environment is encoded in the Lindblad equation eding the localization transition, several numerical ies have found evidence of anomalous subdiffusive sport of particles [18][19][20][21]. The microscopic origin his subdiffusion is still uncertain. The prevailing is that this behavior is caused by Griffiths effects. e, rare regions of strong disorder result in bottlenecks slow transport. Also quasiperiodic system exhibits malous transport, but the correlated nature of the ite potential renders the Griffiths phase explanation an educated guess. The stability of a subdiffusive on, prior to the fully localized phase, is a fascinating sibility, supported by recent numerical works [22,23].
is desirable to find a toy model which capture the key sical aspect of many-body localizing systems, while g treatable. In these notes, we consider the ideas ef. [24] in the non-interacting Aubry-André-Harper H) model. Specifically, the system is evolved coupled to a boundary driving, and to percolated on-site desing. (Each site is subject to dephasing, except for Figure 1. Cartoon summarizing the system when the boundary driving is open. The boundary terms inject/absorb particles from the system. We shall pick i = 0, respectively with probability p, 1 p. The limiting case p = 0 correspond to dephasing at every site, while p = 1 to no dephasing.[To Do: larger fonts]

II. MODELS
The general evolution of open quantum systems coupled to a Markovian environment is encoded in the Lindblad equation    The first term, given by Eq. (6b), describes so-called dephasing processes (energy exchanges with the bath that do not change particle number) and involve the particle density on each site i,n i ≡d † id i . As discussed at the beginning of this section, we take the coupling to the bath γ i to be random and distributed according to Eq. (4).
The last two terms, given in Eq. (6c)-(6d), describe particle injection/ejection contributions acting on the first/last site of the chain, respectively. These processes are necessary to induce a direct current (DC) at late times and therefore to describe transport. We take the coupling to these two baths equal and given by Γ (see Fig. 1).
For brevity, the above characterization is collectively referred to as Markovian Dephasing Bath (MDB). Within this setup, despite the four body interaction due to the dephasing contribution, we can solve the equations of motion for the quantities of interest exactly [53,69] (see Sec. 3.3). We outline the necessary details in App. A.

GOE Bath from Random Regular Graphs
Here we introduce our second setup, which is inspired by the recent proposal of a toy model for the Griffiths phase of the disordered MBL problem [54]. The model is made of M identical copies of Anderson/AAH chains of length L, labeled by the index n = 1, . . . , M . To mimic the effect of thermal inclusions, at each horizontal position i the M sites belonging to different chains are coupled by random hopping terms extracted from a sparse random matrix ensemble, i.e., the ensemble of Random Regular Graphs (RRG) of fixed total connectivity c = 3 [70]. The RRGs are random lattices with a local tree-like structure, loops with typical length of order log M , and no boundary. It is well known that in the absence of disorder the connectivity matrix of a RRG belongs to the Gaussian Orthogonal Ensemble (GOE) [71,72]. Therefore, in the following we refer to this dissipative setting as GOE bath. We present a sketch of this setup in Fig. 2. The Hamiltonian of the model iŝ whered i,n andd † i,n are creation and annihilation operators on the site at position i along the n-th chain, and h i is the inhomogeneous potential (which is identical on all M sites sitting at the same position i). The latter, in accordance with Sec. 3.2.1, we take to be either uniformly distributed in h i ∈ [−W ; W ] or given by a quasi-periodic incommensurate potential, see Eqs. (2) and (3). The parameter t is the hopping rate in the horizontal direction between sites occupying subsequent positions along the chains. γ i is the hopping rate in the transverse direction and it is equal to γ with probability 1 − p and zero with probability p, as in Eq. (4). The notation 〈n, m〉 i indicates pairs of sites n and m with the same horizontal coordinate i connected by an edge of the RRG within the i-th plane.
Note that on each vertical plane a different random realization of the RRG is chosen, in such a way that two sites that are connected by γ within a given layer are (with high probability in the M → ∞ limit) not connected on another layer. This is important as it ensures that the whole lattice can be thought of as an anisotropic random graph, which is locally a tree but has loops whose typical size diverges with the system size. The sites have either connectivity c + 2 if they belong to a layer with γ i = 0 or just 2 in the cases with γ i = 0.
The Anderson insulator with p = 0, corresponding to a uniform dissipative coupling along the chain, has already been studied in Ref. [54]. In this case one finds a localization/delocalization transition when γ becomes of order 1/L. Unlike the MDB set-up described in Sec. 3.2.1 in which each connection to the dephasing bath produces dissipation, in the GOE set-up the RRG couplings are not always effective to ensure dissipation. This implies that while the average number of non-zero γ i couplings along the chain is (1 − p)L, the number of SciPost Phys. 12, 189 (2022) Figure 3: Scheme representing the resonance condition for different layers (labelled i, j, k, l, m, n) with RRG interactions. The density of states in layers i, l and m corresponds to a delta peak, which is characterized by an absence of RRG interactions (γ i = γ l = γ m = 0). On the contrary for layers j, k and n the density of states has a width of 2 c − 1γ as we consider γ j = γ k = γ n = 0. Layers j and k correspond to the scenario in which a resonance is formed, their density of states is non-zero for E = ω and consequently it generates a peak in the LDoS around layers j and k. For layers i, l, m and n however we have no overlap between the probing frequency ω and ρ(E), thus there is no resonance in this case. In general, this means that only the RRG interactions corresponding to the scenario of layers j and k form a resonance in the LDoS. positions in which the GOE-like couplings produces a Local Density of States (LDoS) of order one only scales as ∼ (1 − p eff )L, with p eff ≥ p. This can be understood with simple arguments in the strong disorder limit (W t or λ t depending on the type of on-site potential). When γ = 0 the system corresponds to M identical Anderson insulators with wave-functions strongly localized around the sites i and corresponding eigenvalues close to the local potential −h i . When the RRG couplings are turned on the hopping in the transverse direction lifts up the degeneracy and spreads these M eigenenergies on a semi-circle-like distribution ρ(E) centered around −h i and of width 2 c − 1γ, see App. B. When the probing energy ω falls within this band, a resonance in the LDoS is formed on site i. On the contrary, if ω falls outside this band the RRG couplings do not yield a significant contribution to the LDoS on site i (see Fig. 3 for a sketch). Consequently, the fraction of dephasing planes (1 − p eff ) is roughly given by the probability to have γ i = 0 times the probability that the on-site potential −h i is close enough to the probing energy ω to form a resonance in the LDoS. This means where P(h i ) is the probability distribution of h i and Θ the Heaviside function. More generally, this limiting case (W t or λ t) is a good example to understand why the RRG interactions act like a bath in this set-up. Indeed, as shown in App. B, the presence of these interactions generates a non-vanishing imaginary part in the Green function G(z) = (Ĥ GOE − zˆ ) −1 , with z = ω+i0 + , provided that there is a resonance. This implies that particles on a given site have a finite life-time as they scatter in the system [73], which in other words means that localization is destroyed. With the RRG interactions we thus retrieve a delocalization phenomenon that is usually induced by a bath acting on a system with localized eigenstates.

Observables
In the following we characterize transport and spectral properties of the systems, which are encoded in the frequency-resolved single particle retarded Green's functions whose imaginary part gives the Local Density of States (LDoS). For transport, a natural observable to distinguish between localized, diffusive or intermediate regimes is the resistance of a finite-sized sample, defined as the inverse stationary current flowing through the system, which in the stationary (t → ∞) regime is constant throughout the chain for both protocols (see the App. A).
In the MDB case the average current is defined as The asymptotic current is expected to be independent of i. In the equation above we have introduced the average over the stationary state (infinite time solution of Eq. (5)) 〈•〉 ∞ ≡ tr(•ρ ∞ ).
In the GOE setup instead we directly measure the (zero temperature) dimensionless conductance 1/R GOE at a given energy energy ω through the Fisher-Lee formula [74], which uses the retarded (resp. advanced) Green's function G GOE (ω) (resp. (G GOE (ω)) ) of the system dressed by the leads, and Γ L,R = −2ImΣ L,R as the imaginary part of the self-energies associated with the semi-infinite left and right leads (see App. B for details). Because the GOE case involves the probing energy ω, it is important to note that R MDB and R GOE (ω) are not equivalent: R MDB probes the whole spectrum of the system while R GOE (ω) only characterizes the setup at a given energy. Throughout the rest of the paper we will measure the conductance at ω = 0, around the middle of the band of single-particle eigenstates. Another relevant quantity for our analysis is the single particle retarded Green's function, whose imaginary diagonal part is proportional to the LDoS. Such quantity contains information on the degree of dissipation throughout the system and plays an important role in the theory of localization. We can define the retarded Green's function for the two settings as MDB : respectively. In the MDB setting we have used the time-domain definition of the retarded Green's function, written in terms of the anticommutator of time-evolved fermionic operators and with the theta function ensuring causality. We note that in the context of fermionic open quantum systems the time-evolution has to be performed with the adjoint Lindblad operator as we discuss in App. A.2. In the GOE case it is more convenient to work directly in the frequency representation, z = ω + i0 + , and define the retarded Green's function as the resolvent of the Hamiltonian. As mentioned above, the imaginary part of the retarded Green's function has a direct physical interpretation as LDoS at energy ω, and we study it below. In more explicit form it reads Thanks to the translational invariance within the transverse planes in the GOE setting, the dependence of the Green's functions on the in-layer n indices disappears in the M → ∞ limit. In fact the typical length of the loops in the transverse planes diverges as log M . Therefore, since the local potential is the same on each node of the transverse layer, they all become statistically equivalent and the translational invariance in the transverse direction is recovered. An equivalent expression to the one in Eq. (14) can be obtained for the MDB setting (see e.g. [75]). We note that both the resistance R and the retarded Green's functions G defined above are stochastic variables due to the random coupling to the bath for 0 < p < 1, and have an additional source of randomness for the Anderson model, i.e. the quenched disorder potential h i . We thus study the full probability distribution of these observables. Since the resistance develops a broad distribution with fat tails, P(R) ∼ R −µ , approaching the subdiffusive regime [53], we use its typical value as a proxy characterizing the transport properties. Denoting the average over all sources of disorder as •, we specifically define where we occasionally append the superscript MDB/GOE specifying the protocol under consideration. 2 We also study the exponent µ with which P(R) decays at large R.
We conclude this section with some technical insights on how the calculation of the relevant observables defined above is performed within the two dissipative settings, leaving further details to the Appendices.
In the MDB setting we can compute both the stationary current as well as the retarded Green's function exactly using the equation of motion techniques, as discussed in App. A. This is possible despite the fact that the dephasing jump operators are proportional to the particle density, thus leading to a Lindbladian which is not quadratic in the fermionic operators. The structure of the dissipator however allows for further simplifications: one can indeed show that the equations of motion for single particle correlations, usually coupled to higher order correlators through a full hierarchy, decouple for the dephasing model [53,69,76]. As a result one can obtain closed equations of motion for both single-time and two times single particle correlations, which can be solved in real-space for generic inhomogeneous couplings. Using this approach we can therefore compute both the steady-state current and the LDoS for the MDB, as we review in App. A.
In the RRG/GOE context, we used two approaches. The first one is based on the Cavity Method and has already been explained in Ref. [54], see also App. B. It consists in taking the limit M → ∞ from the start and assuming that loops are so long that they can be neglected. This method is appropriate to compute local observables such as the probability distribution of the LDoS, and leads to recursion relations for the diagonal elements of the Green's functions (given in App. B) which can be easily solved for very long chains.
A second complementary strategy, also explained in App. B, consists in studying the model at finite M coupled to biased reservoirs at its edges and using a relation between the conductivity and the transmission matrix proposed in [74,77]. The coupling to the semi-infinite left and right leads effectively yields a quasi-1d model, which can be solved exactly with the Transfer Matrix method by inverting the full Green's function within each plane by lower-upper (LU) decomposition [78]. This method is appropriate to compute transport properties and global observables such as the conductivity, as it allows one to overcome the drawback of the first approach, which does not account in an exact fashion for all possible paths joining two sites of the first and last external planes which are attached to the leads. However, the exact recursion equations at finite M are much more costly to solve numerically (the time scales as L M 3 ) and we are thus limited to much smaller sizes, L 256, compared to the first strategy. Moreover, keeping M finite leads to extra finite-size effects in the transverse direction. for the statistics of the local density of states. For concreteness, we fix t = 1, Γ = 1/2 and γ = 1/4 for the dephasing protocol, and t = 1, c = 3 and γ = 1/4 for the GOE one.

Resistance, Transport and the Phase diagram
The nature of transport is determined by the scaling of the typical resistance with the system size L. In particular, a power law dependence R typ ∼ L β (16) indicates ballistic transport for β = 0, while diffusion is signaled by β = 1, and corresponds to the Ohm's law [79]. Values of β that are intermediate between 0 and 1 (0 < β < 1) are termed superdiffusive while β > 1 is subdiffusive. The localization limit corresponds to R ∝ exp(L/ξ loc ), with ξ loc the localization length, achieved by a divergent β.
We computed the stationary current as discussed in Sec. 3.3, and we derived the resistance for the two models (see Eqs. (10) and (11)). In Fig. 4 we display the typical resistance R typ , defined in Eq. (15), of the Anderson model with W = 3/2, as a function of system size in double logarithmic scale. We used several values of the probability of decoupling from the local baths, p indicated in the key, in the two protocols, MDB (a) and GOE (b).
In the no-dephasing limit, p = 1, the system exhibits Anderson localization for any W > 0 and the typical resistance grows exponentially with system size. The divergence of β can be observed in the insets where we plot β against p. In the opposite limit, p = 0, in which every site in the chain is subject to dephasing, the system is diffusive R typ ∝ L. This limit is shown with a blue broken line in the main plots and as the limit of the β(p) curves in the inserts. 3 Since the fraction p controls how many sites are unaffected by dephasing, as p becomes large there is an increasing number of insulating intervals. As first discussed in Ref. [23], exponentially distributed rare insulating segments with exponentially large resistance produce subdiffusive transport, R ∼ L β with β > 1, due to Griffiths effects. The trend of the data shown protocols. We obtain the data by extrapolating the system size scaling of the typical resistance R typ (cfr. Eq. (16)). Specifically, we consider system sizes L = 8 ÷ 256 and measure the exponent β = log R typ /log L. The white region correspond to points which are diffusive up to two errorbars. The critical line W c = W (p c ) lies within this area, and signals a transition between a subdiffusive (β > 1) and a diffusive (β = 1) behavior. The black cross markers are examples of extrapolation points, with associated error estimations. The black dashed line p = 1, corresponds to no dephasing. In this limit, the system is always Anderson localized for any W > 0, with a resistance exponentially We note that the GOE bath is less efficient in thermalizing the system compared to the MDB protocol. This is exhibited by the values of the exponent β(p) (insets in Fig. 4): the GOE one reaches anomalous (sub) diffusion values at values of p that are smaller than the ones in the MDB protocol. As detailed in Sec. 3.3, this is due to the fact that the band-width of the RRG perturbation is finite and not all the sites in which γ i > 0 effectively produce dissipation. The value of p in the GOE setting should thus be "renormalized" to p eff ≥ p, as expressed in Eq. (8).
Lastly, we stress that the GOE protocol displays a robust Anderson localization for 1 − p ∼ O(1/L), as schematically depicted in the figure and already discussed in Ref. [54].
From the analysis of the transport properties discussed so far, in particular the evolution of the exponent β in Eq. (16), we can map out a dynamical phase diagram for the Anderson model coupled to the MDB and GOE baths, as a function of p and W , which we plot in Fig. 5 (a)-(b). We see that the crossover from diffusive to sub-diffusive transport extends into a crossover line p c (W ), such that for p < p c (W ) the system is diffusive while for p > p c (W ) it is subdiffusive. The critical dephasing rate decreases with increasing disorder, indicating that a strongly localized Anderson model is more robust to thermal inclusions.
Although, as discussed above, the GOE bath is less efficient to thermalize the samples, the crossover from diffusion to subdiffusion in the two settings is qualitatively very similar.
Finally, we recall that, while the resistance in the MDB setting is computed over all eigenstates, in the GOE case it is computed using only eigenstates around ω = 0. In principle, in the latter case one should also do a finite-size scaling analysis in M . In fact, at any finite M we are considering an effective 1d disordered system which eventually localizes for L M . We can thus expect strong finite size effects and deviations from the power-laws for L M . Such an analysis would demand a heavy numerical work which goes beyond the scope of this paper.

Statistics of the LDoS
We now discuss how the emergence of sub-diffusive transport and Griffiths affects the spectral properties of the system, as encoded in the imaginary part of the retarded local Green's function, i.e. the local Density of States (LDoS).
The statistics of this quantity allows one to capture the dissipation propagation along the chain. In Ref. [54] it was shown that the sub-diffusive phase of the Anderson model with RRG dissipation is characterized by patchworks of insulating segments where ImG i,i (ω) decays exponentially over some length, and rare resonances on which it is of order one. It is therefore interesting to see whether the same physics is captured in the random dephasing Anderson model. For convenience, in the GOE case we fix ω = 0 as a reference energy throughout this subsection. With this choice, in the following we omit the frequency dependence of the LDoS.
In Fig. 6 we plot the histogram of log |Im G i,i | in the Anderson model with MDB (a) and GOE (b) baths. Different curves in the two panels correspond to different values of the parameter p. Concretely, we fix the system size L = 512, and the disorder strength W = 4, and we vary p ∈ [0, 1]. The histograms display a fat tail towards small values of the argument, P( y) ∼ e −τ(p)| y| with τ(p) reported in the inserts of the two graphs. In both cases τ(p) is a decreasing function of p with limit τ(p → 1) → 0. However, the τ(p) takes consistently larger values at p < 1 in the MDB protocol than in the GOE one, due to the fact that in the GOE model the relevant probability of un-coupling to the thermal inclusions is p eff > p. We also note that in the GOE case the histogram shows a bump at large negative values of the argument, which becomes more pronounced as p → 1.
Further insights on the shape of LDoS statistics can be obtained from the analysis of log |Im G i,i | in single samples. More precisely, we now compare the behavior of the Anderson model with the two ways of including dissipation, by using a single and the same realization of {h i } and {γ i }. In Fig. 7 (a) and (b) we display the LDoS along the chain. From these single realizations we can appreciate that for small p (i.e., in the diffusive regime) the LDoS fluctuates around values of order one on all the sites of the chain. Instead, for larger p (i.e. in the In the plots we see, once an effective p eff is considered, the MDB has qualitative features that match those of the GOE protocol. subdiffusive regime) the range of variation of the LDoS is much broader. Figures 7 (a) and (b) demonstrate that the tails of the distribution in Fig. 6, observed at very small values of the LDoS, originate from a very heterogeneous profile of local dissipation in the system. More precisely, we see that very large insulating segments, within which Im G i,i decreases exponentially over the bare localization length (whose size is exponentially distributed) coexist with a few rare resonances where the bath produces dissipation and Im G i,i is of order 1. In App. C we propose a simple model based on this observation to predict the exponent τ(p) characterizing the tail of the LDoS distribution. As shown in the inset of Fig. 6, it agrees well with the numerical results. Furthermore, we can understand the origin of the bump in the very tail of the LDoS distribution for the GOE bath in Fig. 6 as being due to the accumulation of samples without resonances (and thus being fully localized).
Finally, the comparison in Fig. 7 of the LDoS profile between the two types of dissipation confirms what already noted, namely that the MDB bath gives rise to deeper minima in the LDoS profile. The GOE model has larger effective probability of un-coupling to the thermal inclusions, p eff > p. If we take into account this difference and renormalize the coupling γ to the Markovian dissipation, we obtain an overall better qualitative agreement between the two LDoS profiles at a fixed disorder realization as we show in Fig. 7. We conclude that the analysis of the spectral statistics reveals clear signatures of subdiffusive behavior in the Anderson model with the two kinds of dissipation.

The André-Aubry-Harper Chain
In full analogy with the Anderson model analysis, we now investigate the transport and dissipative properties of the AAH model with quasi-periodic disorder and random coupling to a MDB or GOE bath. We recall that in the absence of any coupling to the baths the AAH model displays a non-trivial localization transition at λ c = 2: for λ < λ c the system is ergodic (and  16)) as a function of p. In (c) and (d) λ = 0.4. The asymptotic diffusive limit needs larger systems sizes to establish for increasing p and the long cross-over can be confused with a super-diffusive behavior (see the text for a discussion). The fits are obtained considering the largest system sizes (log L ≥ 4). the transport is ballistic), while for λ > λ c the system is Anderson localized. For p = 1, i.e. in absence of any bulk coupling to dissipation, the model reduces to the one studied in Ref. [80]. Throughout this section we fix t = 1, Γ = 1/2 and γ = 1/4 for the dephasing model, and t = 1, c = 3 and γ = 1/4 for the GOE one.

Resistance, Transport and the Phase Diagram
We start discussing the typical resistance in the steady-state for both the MDB and the GOE AAH models. Our results for λ = 2.4 are given in Fig 8 (a)  First ot all, we see that, in both cases, for sufficiently large system sizes, the typical resistance scales as a power law of system size, i.e. L β for any p < 1, corresponding to a finite coupling to the baths. The exponent β (cfr. Eq. (16)) depends on both the probability p of having a thermal inclusion and the quasi-periodic potential strength λ and it is shown in the inset of Fig. 8.
For λ > λ c , i.e. when the isolated AAH model is in the localized phase, there is a p c = p(λ c ) such that the system is diffusive for 0 ≤ p ≤ p(λ c ), and subdiffusive for p > p(λ c ), see panels (a) and (b) and the insets where the exponent β is seen to deviate from the diffusive limit β = 1 at p c . For λ ≤ λ c , we see that both models show ballistic transport for p → 1, while for p < 1 there is a crossover increasing the system size towards a different scaling behavior which might seem compatible with superdiffusive scaling, L β with β < 1, as recently proposed in Ref. [49,51] for the interacting and isolated AAH. However, being limited to relatively small system sizes, this analysis alone cannot exclude that the apparent superdiffusive behavior observed at intermediate values of L should eventually disappear in the thermodynamic limit, as recently suggested in Ref. [52] for the interacting many-body AAH. We anticipate that the analysis of the LDoS decline any superdiffusion, supporting a diffusive behavior for any p < 1 when λ < λ c (in line with Ref. [16,52]).
Our analysis of the typical resistance is summarized in the transport phase diagram for the AAH model reported in Fig. 9, as a function of the strength of the quasiperiodic potential λ and the probability p of having a dissipative coupling. We find that in both dissipative protocols (MDB on the left and GOE on the right panel) there is a crossover line from subdiffusive to diffusive transport at λ c (p) within the localized phase of the model. On the other hand, as already discussed the transport is diffusive for any p for large enough systems and λ < λ c = 2. Comparing the left panel to the right one we again remark that the GOE bath is less effective in thermalizing the system, with an effective p eff > p.
The emergence of a subdiffusive transport regime in the AAH model under the effect of thermal inclusions, that we found both within the random MDB and the random GOE models, is an intriguing result. To investigate further whether this anomalous transport can be interpreted in terms of Griffiths physics we discuss in the next section the statistics of the LDoS.

Statistics of the LDoS
We conclude the analysis of the AAH model by discussing the statistics of the LDoS across the system evaluated at zero frequency, as done in Sec. 4.2 for the Anderson model. In Fig. 10 we plot the histogram of log |Im G i,i | for the two types of dissipation considered so far and SciPost Phys. 12, 189 (2022)  different values of λ and p. For λ < λ c (panels (a) and (b)) we see that the histogram of the LDoS displays a Gaussian profile, modulo some boundary effects for the MDB, with a variance that increases as p → 1. This provide a robust indication that the system is diffusive, as the initial distribution structure of the h i is washed out by the coupling to the environment. In contrast, this is not the case for λ > λ c (see Fig. 10 (c) and (d)) where a power-law behavior describing the distribution at small values of the LDoS emerges as p is increased. As done for the Anderson Model we can extract an exponent describing the power-law decay at large negative arguments, P( y) ∼ e −τ(p)| y| with τ(p) reported in the inserts of the two graphs. We see that the two exponents are comparable in magnitude, but show an apparent non-linear dependence from p. The presence of a broad distribution for the LDoS is consistent with its behavior in single sample realizations that we had plotted in Fig. 7, showing that for p → 1 both models of thermal inclusions lead for λ > λ c to a very heterogeneous profile of LDoS with few sites having exponentially small dissipation. The overall picture suggests that the AAH model in its localized phase, when coupled to thermal inclusions, also supports a Griffiths like phenomenology, both for what concerns its transport properties, that turn subdiffusive, as well as for the profile of its effective local dissipation.
In the previous sections we discussed transport and dissipation properties in the AAH and Anderson models in the presence of different kinds of thermal inclusions.
One general outcome of this analysis is that the GOE bath is generically less effective in thermalizing the system and restoring diffusion, and leaves a broader signature of anomalous transport effects in its phase diagram, as compared to the MDB setting, as it appears clearly in the transport phase diagrams of Figs. 5-9. We traced back this difference in the spectral properties of the two dissipative settings. In particular, the Markovian dephasing bath, being essentially broadband, is able to establish resonances for any local energy of the system while the RRG/GOE one with its finite bandwidth requires a resonance condition, see Eq. (8).
An interesting aspect of our results are the similarities that emerged in the properties of the Anderson and AAH models coupled to thermal inclusions, which both show a phenomenology compatible with Griffiths effects, in particular for what concerns the subdiffusive transport and the broad statistics of the LDoS. This is an intriguing and somewhat surprising result at first, given the deterministic nature of the quasiperiodic potential which does not obviously lead to rare region effects. One might wonder how much of these similarities has to do with the nature of the coupling to dissipation, which in both Anderson and AAH models we chose to be random. In this respect it is worth emphasizing that the anomalous properties of the AAH model emerged only for λ > λ c , i.e. when in absence of any dissipation the system is in a localized phase, while for λ < λ c the random coupling to thermal inclusions was shown not to be sufficient by itself to destroy diffusion. As such this result seems to point out that the emergence of subdiffusion and broad dissipation statistics is a robust feature of localized phases coupled to thermal inclusions, irrespectively of the nature of the on-site potential.
Finally, it is interesting to comment on the possible connections between our work and a truly interacting disordered many-body problem in the context of MBL. In the random case it was argued in Refs. [53,54] that the role of thermal inclusions could mimic the effect of many-body interactions, which are expected to favor local thermalization. As such the Anderson insulator could be a sensible starting point at strong disorder to understand the destruction of a localized phase due to the proliferation of local ergodic regions. In fact the phenomenology found in both Refs. [53,54] shares many features with the current understanding of the MBL transition. A similar reasoning could be applied to the MBL problem in a quasiperiodic potential. In this respect one could speculate that many-body interactions play a double role, namely to induce local equilibration and at the same time to give rise to a sort of configurational randomness as one could see for example by treating those interactions at the Hartree-Fock level [66,67]. This reasoning justifies our choice of a random coupling to thermal inclusions also in the AAH case. From this perspective one could consider the results discussed in this work also relevant for the MBL/QP problem in a regime where interactions are weak enough to be treated with HF and the disorder is strong enough to be deep in the localized phase and suggest that some sort of subdiffusive scaling could emerge also in that context.

Conclusions
In this work we analyzed the effect of local random dissipation on single particle localized phases of matter, as described by the non-interacting Anderson or AAH models. Several reasons, exposed at the beginning of this manuscript, motivated this study.
In the context of quenched random disorder recent works attempted to give a phenological model of the Griffiths phase found near the MBL transition, as arising from coupling to thermal inclusions. Here, we re-examined two of these works, both of which considered An-derson insulators coupled to different types of dissipative environments, described in terms of random Markovian Dephasing Baths or random GOE baths. In Sec. 4 we provided a thorough comparison of transport and spectral properties of these two settings, thus complementing the existing analysis. As a general outcome we showed that the phenomenology in the two cases is qualitatively similar, with a crossover from sub-diffusive to diffusive transport and the emergence of a broad distribution in the statistics of the local dissipation. Concomitantly, we showed that the GOE bath is typically less effective in thermalizing a localized phase by restoring diffusion, than the Markovian dephasing bath, ultimately due to its finite bandwidth. As a result the properties of an Anderson model coupled to this type of dissipation show anomalous transport and spectral properties in a much broader area of the phase diagram.
A second motivation for this work was to investigate the stability of the localized phase of the quasiperiodic AAH model to thermal inclusions, and to compare its features with the truly random case. In Sec. 5 we extended our analysis to the AAH model coupled to random MDB and GOE dissipation, and we discussed its transport and spectral properties. We provided evidence that even under a quasiperiodic disorder the localized phase is unstable to thermal inclusions towards a phase with Griffiths like phenomenology, including subdiffusive transport and a broad distribution of the local density of states, biased by few rare spots with exponentially small dissipation. Such a phase eventually turns into a conventional diffusive one when the coupling to dissipation becomes more uniform across the chain. On the other hand, we showed that on the delocalized side of the AAH phase diagram this phenomenology is absent and the thermal inclusions lead to diffusive behavior for large enough systems, while finite size samples can display apparent super-diffusive transport. This suggests that the emergence of anomalous transport and spectral properties is a generic feature of localized phases coupled to thermal inclusions, irrespective of the nature of the potential.
In addition to provide two examples of exactly solvable models of disordered systems coupled to random dissipation, both displaying rich physics, our work could be relevant for the MBL problem in the presence of quasiperiodic potentials, as we discuss in Sec. 6. In this respect, it would be interesting to compare more closely the phenomenology that emerges from our toy models with results obtained for the fully interacting case.

A.1 Stationary correlation functions and current
To evaluate the current we only need the evolution of the (instantaneous) spatial correlation matrix In the above equation, we denote time t. By differentiating with respect to the time variable t and using the Lindblad equation (5) we find Remarkably, despite the fact that dephasing acts as a four body interaction, the equations of motion for the two-point functions are closed and can be solved exactly. We also note that these equations are decoupled from the particle-pair creation/annihilation correlations. In Eq. (18) we implicitly imposed open boundary condition C 0,i (t) = C i,L+1 (t) = 0 for all i. The stationary state is obtained by solving ∂ t C i, j = 0 for all i, j. For this purpose, it is convenient to rephrase Eq. (18) in matrix form. We introduce the matrix capturing the dissipation contributions with elements D i,k = δ i,k (γ i +δ i,1 Γ +δ i,L Γ ), and the matrix = ( ), with components P i,k (t) = 2δ i,k (γ k C k,k (t) + Γ δ i,L δ k,L ) ≡ p k (t)δ i,k . With these definitions, we recast (18) as where we avoided writing the time dependencies. The infinite time solution is obtained with standard ordinary differential equation tools [53], and reads We note that this is a self-consistency equation due to the dependency of on . This expression can be further massaged using the spectral decomposition of the matrix = p λ p |φ R (p)〉〈φ L (p)| , 〈φ L (q)|φ R (p)〉 = δ p,q , where with a slight abuse of notation we used the Dirac notation for the L-dimensional complex vectors |φ R/L (p)〉 = i φ which we use in Eq. (20) to derive the self-consistent equation [53,80] C i, j (∞) = 2Γ Θ i, j,L + 2 k γ k Θ i, j,k C k,k (∞) , Lastly, the stationary current is j ∞ ≡ j i (∞) = 2Im C i−1,i (∞). We conclude this subsection by remarking that not all the components of the tensor Θ are needed to extract the stationary value of the current, which is constant throughout the system. SciPost Phys. 12, 189 (2022) where G (r) i is an M × M complex symmetric matrix (and similarly for the measure in absence of the left neighboring layer). By Gaussian integration, it is straightforward to show that the recursion relations for the left and right elements of the measure read: Here, C (i) nm is the connectivity matrix of the RRG on the i-th layer whose elements are equal to 1 if n and m are connected by the RRG and zero otherwise, and γ i = γ with probability 1 − p and zero with probability p. Considering an open system the boundary conditions on i = 1 and i = L are The first two equations for the Green's functions can be easily solved by inverting the matrices on the right hand side by lower-upper (LU) decomposition [78] on each layer. Note that these equations are exact but can only be solved for finite M .
However, in the M → ∞ limit a drastic simplification arises. In this limit, the typical length of the loops in the transverse planes diverges and the Gaussian measure on the c neighbors of a given site within a given layer factorizes in absence of the central site, since these c neighbors belong to disconnected branches of the graph and are uncorrelated. Moreover, as explained in the main text, since the local potential is the same on all the vertices of the transverse planes, in the M → ∞ limit they all become statistically equivalent and a translational invariance within each plane is restored. Equations (31) and (32) are thus drastically simplified. Following Ref. [54], one can show that: with c the uniform connectivity of the RRG on the planes, and where G (r,l,v) i are the diagonal elements of the Green's function on any site of the i-th layer in absence of its right neighbor, its left neighbor, and one of its c neighbors in the transverse plane, respectively. Once these so-called "cavity" Green's functions are known on each layer, one can finally compute the diagonal elements of the Green's function of the original problem on layer i (see Ref. [54] for more details): via We insist upon the fact that there is no dependence on the in-layer indices n, m within this approximation. Besides, focusing on the limiting case of independent layers (t = 0) we can exactly obtain the density of states in each layer. As a matter of fact for this case we can solve Eq. (35) and obtain around e −l/2ξ loc (the minimum value of the density being attained roughly at the middle of the segment) to O(1). Following the analysis in [53] (in particular the third section) the total number of layers with log |Im(G i,i )| equal to −α (up to sub-leading terms) is where n(l) is the number of isolating segments of length l in the system. Equation (43) simply translates the fact that only isolating segments longer than 2ξ loc α contain layers where log |Im(G i,i )| = −α: these layers being approximately two in number, one when the LDoS decreases exponentially with i, and one when the LDoS increases exponentially with i. Then, the probability distribution is given by In fact, if we consider that each layer in the system has a probability p to form a local resonance in the LDoS (ImG i,i = O(1)) it can be proven that n(l) is a random variable following the Poisson distribution where 〈•〉 indicates the average over all variables {n(l)} l∈ . For p close to one this simple approach predicts a tail τ(p) in the probability distribution 〈P log |Im(G i,i )| = −α]〉 ∼ e −ατ(p) , with τ(p) = −2ξ loc log(p) .
When compared to our numerical results we observe a good agreement with this prediction, see the inset in Fig. 6.