Dephasing-induced mobility edges in quasicrystals

In this section we derive the classical master equation and explicit form of the Markov transition matrix $W$, as given by Eqs.(4) and (5) of the main text, in the limit of a large dephasing rate.
In the single-particle sector of Hilbert space, the Lindblad master equation yields the following evolution equations for the density matrix elements $\rho_{n,m}$ ME1 ; ME2

$$\frac{d\rho_{n,m}}{dt^{\prime}}=(i/\gamma)(J_{n-1}\rho_{n-1,m}+J_{n}\rho_{n+1,m})-(i/\gamma)(J_{m}\rho_{n,m+1}+J_{m-1}\rho_{n,m-1})+(i/\gamma)(V_{m}-V_{n})\rho_{n,m}-(1-\delta_{n,m})\rho_{n,m}$$

where $t^{\prime}=\gamma t$ is the dimensionless time scaled to the coherence time $1/\gamma$. When the decay rate $\gamma$ of coherences is much larger than any terms $|J_{n}|$ and $|V_{n}|$, the quantities $|J_{n}|/\gamma$ and $|V_{n}|/\gamma$ on the right hand side of the above equation are small, and we highlight this fact by writing the equation in the form

$$\frac{d\rho_{n,m}}{dt^{\prime}}=\epsilon(i/\gamma)(J_{n-1}\rho_{n-1,m}+J_{n}\rho_{n+1,m})-\epsilon(i/\gamma)(J_{m}\rho_{n,m+1}+J_{m-1}\rho_{n,m-1})+\epsilon(i/\gamma)(V_{m}-V_{n})\rho_{n,m}-(1-\delta_{n,m})\rho_{n,m}$$

(S1)

where $\epsilon$ is a flag that marks the order of magnitude of such small quantities. In such form, one can solve perturbatively the equation by a standard multiple-time scale method (see for instance multiple ), and at the end of the calculation we can set $\epsilon=1$. After introduction of the multiple time scales $T_{0}=t^{\prime}$, $T_{1}=\epsilon t^{\prime}$, $T_{2}=\epsilon^{2}t^{\prime}$, … we can search for a solution of Eq.(S1) as a power series expansion in $\epsilon$, i.e.

$$\rho_{n,m}(t^{\prime})=\rho_{n,m}^{(0)}(t^{\prime})+\epsilon\rho_{n,m}^{(1)}(t^{\prime})+\epsilon^{2}\rho_{n,m}^{(2)}(t^{\prime})+...$$

(S2)

The power series expansion is meaningful provided that each term in the expansion does not show secularly growing terms in time, which would prevent the validity of Eq.(S2) after a short evolution time. In order to avoid the appearance of secular terms, solvability conditions should be satisfied, which will provide the dynamical evolution of the leading-order density matrix elements $\rho_{n,m}^{(0)}$ on the different time scales $T_{0}$, $T_{1}$, $T_{2}$, …. Substitution of Eq.(S2) into Eq.(S1) and using the derivative rule $d/dt^{\prime}=\partial_{T_{0}}+\epsilon\partial_{T_{1}}+\epsilon^{2}\partial_{T_{2}}+...$, after collecting the terms of the same order in $\epsilon$ a hierarchy of equations for successive corrections to $\rho_{n,m}(t^{\prime})$ is obtained. At leading order $\sim\epsilon^{0}$ one has

$$\frac{\partial\rho_{n,m}^{(0)}}{\partial T_{0}}=-(1-\delta_{n,m})\rho_{n,m}^{(0)}.$$

(S3)

For $n=m$, the solution to Eq.(S3), which does not display secular term growing as $\sim T_{0}$, is simply given by $\rho_{n,n}^{(0)}=P_{n}(T_{1},T_{2},...)$, i.e. it is independent of $T_{0}$. This means that the population dynamics evolves on the slow time scales $T_{1}$, $T_{2}$, …, i.e. the relaxation toward the stationary maximally-mixed state of the Liouvillian occurs on such slow time scales. For $n\neq m$, the solution to Eq.(S3) is given by $\rho_{n,m}^{(0)}=\rho_{n,m}^{(0)}(0)\exp(-T_{0})$, indicating that any initial coherence in the system is rapidly damped on the fast time scale $T_{0}$. Since the relaxation dynamics of populations occurs on the slow time scales, we can thus disregard such a fast transient decay of coherences, and let at leading order

$$\rho_{n,n}^{(0)}=P_{n}(T_{1},T_{2},...)\;,\;\;\rho_{n,m}^{(0)}=0\;\;(n\neq m).$$

(S4)

The change of populations $P_{n}$ occurs on the longer time scales $T_{1}$, $T_{2}$, … due to incoherent hopping mediated by the small values of coherences at higher orders, and can be calculated from the solvability conditions by pushing the analysis up to order $\sim\epsilon^{2}$.
At order $\sim\epsilon$ one obtains

$$\partial_{T_{0}}\rho_{n,m}^{(1)}+(1-\delta_{n,m})\rho_{n,m}^{(1)}=G_{n,m}^{(1)}$$

(S5)

where we have set

$$G_{n,m}^{(1)}\equiv(i/\gamma)(J_{n-1}\rho_{n-1,m}^{(0)}+J_{n}\rho_{n+1,m}^{(0)})-(i/\gamma)(J_{m}\rho_{n,m+1}^{(0)}+J_{m-1}\rho_{n,m-1}^{(0)})+(i/\gamma)(V_{m}-V_{n})\rho_{n,m}^{(0)}-\frac{\partial\rho_{n,m}^{(0)}}{\partial T_{1}}$$

(S6)

The solution to Eq.(S5) reads

	$\displaystyle\rho_{n,n+1}^{(1)}$	$\displaystyle=$	$\displaystyle\frac{iJ_{n}}{\gamma}(P_{n+1}-P_{n})$
	$\displaystyle\rho_{n,n-1}^{(1)}$	$\displaystyle=$	$\displaystyle-\frac{iJ_{n-1}}{\gamma}(P_{n}-P_{n-1})$		(S7)
	$\displaystyle\rho_{n,m}^{(1)}$	$\displaystyle=$	$\displaystyle 0\;\;m\neq n\pm 1$

with the solvability conditon

$$\frac{\partial P_{n}}{\partial T_{1}}=0.$$

(S8)

Equation (S7) provides the leading-order terms ($\sim\epsilon$) of non-vanishing coherences in the system, which occur along the two diagonals $m=n\pm 1$.
At order $\sim\epsilon^{2}$ one finally obtains

$$\partial_{T_{0}}\rho_{n,m}^{(2)}+(1-\delta_{n,m})\rho_{n,m}^{(2)}=G_{n,m}^{(2)}$$

(S9)

where we have set

$$G_{n,m}^{(2)}\equiv(i/\gamma)(J_{n-1}\rho_{n-1,m}^{(1)}+J_{n}\rho_{n+1,m}^{(1)})-(i/\gamma)(J_{m}\rho_{n,m+1}^{(1)}+J_{m-1}\rho_{n,m-1}^{(1)})+(i/\gamma)(V_{m}-V_{n})\rho_{n,m}^{(1)}-\frac{\partial\rho_{n,m}^{(0)}}{\partial T_{2}}$$

(S10)

In particular, if we specialize Eq.(S9) for the populations ($m=n$), one obtains the solvability condition $G_{n,n}^{(2)}=0$, i.e.

$$\frac{\partial P_{n}}{\partial T_{2}}=(i/\gamma)(J_{n-1}\rho_{n-1,n}^{(1)}+J_{n}\rho_{n+1,n}^{(1)})-(i/\gamma)(J_{n}\rho_{n,n+1}^{(1)}+J_{n-1}\rho_{n,n-1}^{(1)}).$$

(S11)

Substitution of Eq.(S7) into Eq.(S11) yields

$$\frac{\partial P_{n}}{\partial T_{2}}=-\frac{2}{\gamma^{2}}(J_{n}^{2}+J_{n-1}^{2})P_{n}+\frac{2J_{n}^{2}}{\gamma^{2}}P_{n+1}+\frac{2J_{n-1}^{2}}{\gamma^{2}}P_{n-1}.$$

(S12)

Once the solvability condition Eq.(S12) is satisfied, the higher-order corrections to coherences, at order $\sim\epsilon^{2}$, can be obtained by solving Eq.(S9).
If we stop the analysis at this order $\sim\epsilon^{2}$, using the derivative rule $dP_{n}/dt=\gamma(\partial_{T_{0}}P_{n}+\epsilon\partial_{T_{1}}P_{n}+\epsilon^{2}\partial_{T_{2}}P_{n})=\gamma\epsilon^{2}\partial_{T_{2}}P_{n}$ in physical time $t$, after letting $\epsilon=1$ from Eq.(S12) one then obtains Eqs.(4) and (5) given in the main text.

Let us consider the single-particle coherent evolution in the tight-binding lattice and let us assume that at every time interval $\Delta t$ the phase of the wave function $\psi_{n}$ at any lattice site $n$ is randomized. Such a randomized phase process basically emulates dephasing effects in the dynamics, and provides a simple tool to experimentally emulate dephasing phenomena in photonic quantum walks S1 . The time evolution of the wave function amplitudes $\psi_{n}(t)$ read

$$i\frac{d\psi_{n}}{dt}=\sum_{m}H_{n,m}\psi_{m}+\psi_{n}\sum_{\alpha=1,2,3,...}\varphi_{n}^{(\alpha)}\delta(t-\alpha\Delta t)$$

(S13)

where $H_{n,m}$ are the matrix elements of the single-particle tight-binding Hamiltonian $H$ and the last term on the right hand sides of Eq.(S13) describes the dephasing process. Here $\varphi_{n}^{(\alpha)}$ are assumed to be uncorrelated stochastic phases, both in site index $n$ and time step $\alpha$, with a given probability density function. Fully coherent dynamics is obtained by letting $\varphi_{n}^{(\alpha)}=0$, whereas fully incoherent (classical) dynamics is obtained by assuming a uniform distribution in the range $(-\pi,\pi)$ for the probability density function. Under fully coherent dynamics, the wave function amplitudes evolve according to $\psi_{n}(t)=\sum_{m}U_{n,m}(t)\psi_{m}(0)$, where $U(t)=\exp(-iHt)$ is the coherent propagator. On the other hand, for fully incoherent dynamics indicating by $P_{n}(t_{\alpha})={\overline{|\psi_{n}(t_{\alpha})|^{2}}}$ the occupation probabilities (populations) at various sites of the lattices, where $t_{\alpha}=\alpha\Delta t$ and the overbar denotes statistical average, the time evolution is described by the classical map

$$P_{n}(t_{\alpha+1})=\sum_{m}\mathcal{U}_{n,m}P_{m}(t_{\alpha})$$

(S14)

where $\mathcal{U}_{n,m}=|U_{n,m}(\Delta t)|^{2}$ in the incoherent propagator. Equation (S14) can be readily obtained by solving Eq.(S13) in each time interval $\Delta t$ and then taking the statistical average. The classical master equation (4) given in the main text is finally obtained in the small $\Delta t$ limit. In fact, assuming a short time interval $\Delta t$ between successive stochastic phases, one can expand the coherent propagator in power series of $\Delta t$, yielding $U(\Delta t)\simeq 1-i\Delta tH-(1/2)\Delta t^{2}H^{2}$ and thus

$$\mathcal{U}_{n,m}=\delta_{n,m}+\Delta t^{2}\left(|H_{n,m}|^{2}-\delta_{n,m}\sum_{l}|H_{n,l}|^{2}\right)$$

(S15)

up to the order $\sim\Delta t^{2}$. Since $P_{n}(t_{\alpha})$ evolves slowly with index $\alpha$, i.e. after each short time interval $\Delta t$, one can set $P_{n}(t+\Delta t)=P_{n}(t)+(dP_{n}/dt)$ and consider $t_{\alpha}=t$ a continuous variable. From Eqs.(S14) and (S15) one then obtains the classical master equation

$$\frac{dP_{n}}{dt}=\sum_{m=1}^{N}W_{n,m}P_{m}(t)$$

(S16)

where the elements of the Markov transition matrix $W$ are given by

$$W_{n,m}=\Delta t\left(|H_{n,m}|^{2}-\delta_{n,m}\sum_{l}|H_{n,l}|^{2}\right)=\left\{\begin{array}[]{cc}-\Delta t\sum_{l\neq n}|H_{n,l}|^{2}&n=m\\ \Delta t|H_{n,m}|^{2}&n\neq m\end{array}\right.$$

(S17)

For the tight-binding model describing the quasicrystal, one has $H_{n,m}=V_{n}\delta_{n,m}+J_{n}\delta_{n,m-1}+J_{n-1}\delta_{n,m+1}$, so that one obtains

$$W_{n,m}=\left\{\begin{array}[]{cc}-\Delta t(J_{n}^{2}+J_{n-1}^{2})&n=m\\ \Delta tJ_{n}^{2}&n=m-1\\ \Delta tJ_{n-1}&n=m+1\\ 0&n\neq m,m\pm 1\end{array}\right.$$

(S18)

which reduces to Eq.(5) of the main text provided that we assume $\Delta t=2/\gamma$.

The localization properties of a wave function $\psi_{n}^{(l)}$ with energy $E_{l}$, normalized as $\sum_{n=1}^{L}|\psi_{n}^{(l)}|^{2}=1$, are characterized by the generalized inverse participation ratio, defined by

$${\rm IPR}_{l}^{(q)}\sum_{n=1}^{L}|\psi_{n}^{(l)}|^{2q},$$

(S19)

where $q>0$ is a positive number, and by the exponent $\beta_{l}^{(q)}$, which is given by (see for instance r1 ; r2 ; r3 ; r4 ; r5 )

$$\beta_{l}^{(q)}=\lim_{L\rightarrow\infty}\frac{\ln{\rm IPR}_{l}^{(q)}}{\ln(1/L)}.$$

(S20)

The quantity $D_{l}(q)=\beta_{i}^{(q)}/(q-1)$ is known as the fractal dimension. Note that for $q=2$ one has $D_{l}(q)=\beta_{l}^{(q)}$. For extended (ergodic) and localized states one has $D_{l}(q)=1$ and $D_{l}(q)=0$, respectively, whereas any other behavior of $D_{l}(q)$ implies multifractality r5 , i.e. critical states. In the finite-size scale analysis of localization, the golden ratio $\alpha=(\sqrt{5}-1)/2$ is approached by the Fibonacci numbers via the relation $\alpha=\lim_{l\rightarrow\infty}F_{l-1}/F_{l}$, where the Fibonacci numbers $F_{l}=1,1,2,3,5,8,13,21,34,55,89,144,233,377,$ $610,987,1597,2584,4181,6765,10946,17711,...$ are defined recursively by $F_{l+1}=F_{l-1}+F_{l}$, with $F_{0}=F_{1}=1$. Numerically, we successively change the system size $L=F_{l}$ to approach the irrational number, and extrapolate IPR${}_{l}^{(q)}$ and $\beta_{l}^{(q)}$ by their asymptotic behavior as $L$ is increased. Open boundary conditions are assumed in the simulations, with system size $L$ typically spanning the range of Fibonacci numbers from $L=34$ to $L=17711$.

In the off-diagonal Aubry-André model with dephasing, the Markov transition matrix $W$ displays incommensurate disorder in both hopping rates and on-site energies, and its energy spectrum shows a typical Cantor-like set structure of incommensurate models r6 , with infinitely many small bands separated by small gaps. The behavior of energy spectrum of $W$ versus $\kappa=B/A$, computed by diagonalization of the Markov matrix $W$ on a lattice of size $L=4181$, is shown in Fig.S1. The spectrum basically comprises three main pseudo bands, labelled as I,II and III in the figure, separated by large gaps. For $\kappa$ below the critical value $\kappa_{c}\simeq 0.4$, all eigenvectors of $W$ are extended, whereas for $\kappa>\kappa_{c}$ coexistence of extended and localized eigenstates in suggested from the IPR behavior, which is indicated by the different colors in the plot.
To get deeper insights into the localization properties of the eigenvectors and the localization-delocalization transition, extended numerical simulations have been performed computing the behavior of Lyapunov exponent (inverse of localization length), fractal dimension and integrated level spacing distribution (ILSD). Typical numerical results are shown in Figs.S2 and S3 for $\kappa$ well above ($\kappa=0.5$, Fig.S2) and just above ($\kappa=0.407$, Fig.S3) the critical point $\kappa_{c}\simeq 0.4$. In the former case (Fig.S2) all eigenvectors in pseudo bands I and II are localized, while in pseudo band III they are all extended, as clearly shown in Figs.S2(b) and (c). The localization of the eigenvectors is characterized by the Lyapunov exponent (inverse of the localization length), which is depicted in Fig.S2(c). The Lyapunov exponent $\gamma(\lambda_{l})$ for the eigenvector $\psi_{n}^{(l)}$ of $W$ with eigenvalue (energy) $\lambda_{l}$ has been numerically computed using the relation r7 ; r8 ; r9

$$\gamma(\lambda_{l})=\lim_{L\rightarrow\infty}\frac{1}{L-1}\sum_{\lambda_{n}\neq\lambda_{l}}\log|\lambda_{n}-\lambda_{l}|-\lim_{L\rightarrow\infty}\frac{1}{L-1}\sum_{n=1}^{L-1}\log|W_{n,n+1}|$$

(S21)

assuming a large lattice size (typically $L=6765$). For an energy $\lambda$ in the spectrum, $\gamma(\lambda)$ vanishes for an extended or critical state, whereas it assumes a finite value for an exponentially-localized state, $1/\gamma(\lambda)$ providing the characteristic localization (decay) length of the eigenvector. As shown in Fig.S2(c), all eigenvectors in spectral region I display a finite Lyapunov exponent, with a localization length that increases as $\lambda$ is increased but that remains finite over the entire spectral region. Conversely, in the pseudo bands II and III the Lyapunov exponent is very close to zero, indicating that in such regions all the eigenvectors are either extended or critical. Since the Lyapunov exponent jumps from finite to almost vanishing values, the mobility edge for $\kappa=0.5$ falls in the wide gap separating the pseudo bands I and II. This behavior also suggests that there are not critical states, and the mobility edge separates extended (ergodic) from exponentially-localized eigenvectors. A few examples of localized and extended eigenstates, corresponding to the four eigenvalues $\lambda_{1}=-0.093$, $\lambda_{2}=-0.0329$, $\lambda_{3}=-4\times 10^{-3}$ and $\lambda_{4}=-2\times 10^{-4}$, are depicted in Fig.S2(d). To check the absence of critical states, Fig.S2(e) shows the behavior of ${\rm ln}$(IPR) (with $q=2$) versus lattice size ${\rm ln}(1/L)$ for the four eigenvectors of panel (d). The estimated slopes $\beta$ of the four curves, $\beta\simeq 0$ for the eigenvector with eigenvalue $\lambda=\lambda_{1}$ and $\beta\simeq 1$ for the eigenvectors with eigenvalues $\lambda=\lambda_{2},\lambda_{3}$ and $\lambda_{4}$, indicate that the eigenvectors are either localized ($\lambda=\lambda_{1}$) or extended ergodic ($\lambda=\lambda_{2},\lambda_{3},\lambda_{4}$), i.e. they do not display multifractality. The nature of the eigenstates is also confirmed by the behavior of exponents $\beta^{(q)}$ versus $q$, for a fixed value of lattice size $L$, which is depicted in Fig.S2(f).

A different behavior is found near the critical point (Fig.S3). In this case, besides extended and localized eigenstates, there is evidence of critical states at the mobility edge separating extended and localized states, which falls in a narrow sub-band of pseudo band I [inset of Fig.S3(a)]. The eigenvector with eigenvalue $\lambda=\lambda_{3}$, depicted in the third row of Fig.S3(d), is neither fully localized nor fully extended, i.e. it is a critical state displaying multifractality, as indicated by the behavior of $\beta$ for this eigenvector shown in Figs.S3(e) and (f).

Finally, we performed an analysis of level spacing statistics in order to check the appearance of critical states. Level spacing statistics is an important tool in the study of spectra of disordered systems, which can provide major insights into the localization properties and critical points of a system r10 ; r11 . It is well known that for a quantum system with uncorrelated disorder localized and delocalized levels follow different universal level spacing distributions $P(s)$ in the thermodynamic limit, defined as the probability density of level spacings $s$ of the adjacent levels r10 ; r11 . In the extended phase $P(s)$ follows the Wigner surmise $P_{W}(s)\sim s\exp(-cs^{2})$ distribution (with $c$ a constant), whereas in the localized phase $P(s)$ follows a Poisson law $P_{P}(s)\sim\exp(-s)$. Thus, as disorder increases, the metal-insulating transition is accompanied by a transition of the level spacing distributions from the Wigner surmise to the Poisson law. The Wigner surmise is characterized by the typical level repulsion since $P(s)\rightarrow 0$ as $s\rightarrow 0$, whereas for the Poisson distribution level repulsion is absent. In a disordered system displaying mobility edges, the level spacing distribution $P(s)$ changes discontinuously at the mobility edge from $P_{P}(s)$ to $P_{W}(s)$. In models with incommensurate potentials, the most notable one being the Aubry-André model, the energy level spacing distributions in the extended and localized phases do not follow such universal behaviors r12 ; r13 ; r14 ; r15 ; r16 . Interestingly, at the critical point, where the eigenvectors are critical, $P(s)$ is singular as $s\rightarrow 0$ and displays the inverse power law $P(s)\sim s^{-3/2}$, indicating energy level clustering (rather than repulsion). Such an inverse power law and level clustering have been suggested as signatures of critical wave functions in incommensurate potential models r13 ; r14 ; r15 . Owing to the fragmented structure of energy spectrum, formed by an uncountable set of levels (like in the Cantor set), rather than computing the level spacing distributions $P(s)$ it is more convenient to calculate an integrated level-spacing distribution (ILSD) which, apart from normalization reads r13 ; r14 ; r15

$${\rm ILSD}(s)=\int_{s}^{\infty}ds^{\prime}P(s^{\prime})$$

(S22)

and provides the fraction of the total number of gaps larger than some size $s$. The ILSD can be computed in any given range of energies, either in the localized or extended regions, to characterize the localization features of the eigenvectors. For a given range of energies $(\mathcal{E}_{1},\mathcal{E}_{2})$, the level spacing $s$ is defined as $s=(E_{n+1}-E_{n})/(W/N_{E})$, where $E_{1}\leq E_{2},...,E_{N_{E}}$ are the $N_{E}$ energy levels that fall in the range $(\mathcal{E}_{1},\mathcal{E}_{2})$ and $W=(\mathcal{E}_{2}-\mathcal{E}_{1})$. The ILSD can be normalized by introducing a lower cutoff $s_{0}>0$ r13 . In Figs.S2(g) and S3(g) we show the numerically-computed behavior of the ILSD in the spectral intervals I,II and III defined by the three main pseudo bands. A clear different behavior is observed for the ILSD curves in region I for $\kappa=0.5$, where the ME falls in a gap and there are not critical states [Fig.S2(g)], and $\kappa=0.407$, where the ME falls in a narrow sub-band of region I [Figs.S3(g)]. In the latter case the ILSD curve displays a steep increase near $s=0$ [see inset in Fig.S3(g)], indicating level clustering. Such a result is a signature of the emergence of critical states for $\kappa=0.407$, which are absent when $\kappa=0.5$.

The key result of dephasing-induced ME presented in the main manuscript for the off-diagonal Aubry-André model has been demonstrated in the strong dephasing regime, which justifies a classical reduction of the quantum master equation Kamia . Here we show that such a phenomenon persists for finite values of the dephasing rate as well, i.e. beyond the classical limit. This entails the diagonalization of the full Liouvillian superoperator $\mathcal{L}$ entering in the Lindblad master equation. In a lattice of size $L$ the density matrix $\rho$ in the single particle sector of Hilbert space is represented by a vector of size $L^{2}$ and the Liouvillian superoperator $\mathcal{L}$ by a $L^{2}\times L^{2}$ square matrix. This makes the diagonalization problem computationally challenging for very large system sizes, however one can provide some insightful results of the cross-over between fully coherent and fully incoherent regimes taking a relatively large value of $L$ to emulate a quasicrystal, yet keeping the eigenvalue problem tractable with standard MatLab eig solver. In the numerical simulations, the golden mean $\alpha$ has been approximated by the Fibonacci ratio $\alpha=34/55$, corresponding to a system size $L=55$. The eigenvalue problem reads $\mathcal{L}\rho=\lambda\rho$, and the stationary state $\rho_{s}$, corresponding to the vanishing eigenvalue $\lambda=0$, is the maximally-mixed states, with vanishing coherences and populations uniformly distributed over the $L$ sites of the lattice. The other eigenvectors of the Liouvillian superoperator correspond to eigenvalues $\lambda$ with negative real part, and only those with long lifetime (i.e. small values of $|{\rm Re}(\lambda)|$) are relevant for the relaxation dynamics. For a given eigenvector $\rho$ of $\mathcal{L}$, the IPR is defined for the populations as IPR$=\sum_{n}\rho_{n,n}^{2}$. An example of the spectrum of the Liouvillian superoperator $\mathcal{L}$ and corresponding behavior of IPR of all eigenvectors, for parameter values $A=1$, $B=0.7$ and $\gamma=A=1$, is given in Fig.S4. The eigenvalue spectrum of $\mathcal{L}$ comprises a majority of eigenvalues with real part grouped near $-\gamma$, corresponding to the most damped eigenvectors containing the coherences, and a set of eigenvalues that condensate near $\lambda=0$, i.e. closer to the stationary maximally-mixed eigenvector $\rho_{s}$ of $\mathcal{L}$. The behavior of IPR clearly indicates that for such set of eigenvectors the populations are markedly localized in the region below the mobility edge $\lambda_{m}$, and extended above such a ME. Hence we have clear indications that ME appear even for weak-to-moderate values of the dephasing rate $\gamma$. To corroborate such an indication, in Fig.S5 we plot the behavior of IPR_min and IPR_max (among all eigenvectors of $\mathcal{L}$) versus $\kappa=B/A$ for a few increasing values of the dephasing rate $\gamma$. The figure clearly shows that there is still a phase transition at some critical value $\kappa_{c}$ even for weak values of the dephasing rate, with $\kappa_{c}$ approaching the value $\kappa_{c}=0.4$ in the strong dephasing (classical) regime $\gamma/A\gg 1$ and the value $\kappa_{c}=1$ in the weak dephasing regime $\gamma/A\ll 1$, as expected.

Here we show that the off-diagonal Aubry-André model can be emulated by the discrete-time quantum walk of photons in a synthetic mesh lattice realized by two coupled fiber loops with slightly unbalanced lengths. The light pulse dynamics with inhomogeneous coupling angles $\beta_{n}$ is described by the set of discrete-time equations S2

	$\displaystyle u^{(m+1)}_{n}$	$\displaystyle=$	$\displaystyle\left(\cos\beta_{n+1}u^{(m)}_{n+1}+i\sin\beta_{n+1}v^{(m)}_{n+1}\right)\exp(i\phi_{n}^{(m)})\;\;\;\;\;\;$		(S23)
	$\displaystyle v^{(m+1)}_{n}$	$\displaystyle=$	$\displaystyle i\sin\beta_{n-1}u^{(m)}_{n-1}+\cos\beta_{n-1}v^{(m)}_{n-1}$		(S24)

where $u_{n}^{(m)}$ and $v_{n}^{(m)}$ are the pulse amplitudes at discrete time step $m$ and lattice site $n$ in the two fiber loops, $\beta_{n}$ is the site-dependent coupling angle, and $\phi_{n}^{(m)}$ are uncorrelated stochastic phases with uniform distribution in the range $(-\pi,\pi)$.

The coherent quantum walk is obtained by Eqs.(S23) and (S24) by letting $\phi_{n}^{(m)}=0$. Let us assume a coupling angle $\beta_{n}$ close to $\pi/2$, i.e. let us assume

$$\beta_{n}=\frac{\pi}{2}-\theta_{n}$$

(S25)

with $|\theta_{n}|\sim\epsilon\ll 1$. Up to order $\sim\epsilon$, we may set $\sin\beta_{n}\simeq 1$ and $\cos\beta_{n}\simeq\theta_{n}$ in Eqs.(S23) and (S24). After letting $u_{n}^{(m)}=(i)^{m}F_{n}^{(m)}$, elimination of the amplitudes $v_{n}^{(m)}$ from Eqs.(S23) and (S24) yields the following second-order difference equation for $F_{n}^{(m)}$

$$F^{(m+1)}_{n}=F^{(m-1)}_{n}-i\left(\theta_{n}F_{n-1}^{(m)}+\theta_{n+1}F_{n+1}^{(m)}\right)$$

(S26)

which is accurate up to order $\sim\epsilon$. Equation (S26) indicates that the amplitude $F_{n}^{(m)}$ changes by a small quantity, of order $\epsilon$, every two discrete time steps, i.e. when $m$ is increased to $(m+2)$. Therefore, the most general solution to Eq.(S26) can be written as the superposition

$$F_{n}^{(m)}=\psi_{n}^{(+)}(m)+(-1)^{m}\psi_{n}^{(-)}(m)$$

(S27)

where the amplitudes $\psi_{n}^{(\pm)}(m)$ are slowly-varying functions of $m$ and satisfy the decoupled equations

$$i\frac{\partial\psi_{n}^{(\pm)}}{\partial m}=\pm\left(J_{n}\psi_{n+1}^{(\pm)}+J_{n-1}\psi_{n-1}^{(\pm)}\right)$$

(S28)

where we have set

$$J_{n}=\frac{1}{2}\theta_{n+1}\simeq\frac{1}{2}\cos\beta_{n+1}$$

(S29)

and treated $m$ as continuous variable. Clearly, Eq.(S28) corresponds to the single-particle Schrödinger equation in a one-dimensional lattice with inhomogeneous hopping rates $J_{n}$. The off-diagonal Aubry-André model is thus retrieved by letting $\beta_{n}=\pi/2-2A-2B\cos(2\pi\alpha n)$, with $A,B\ll 1$.