Localization of light in three dimensions: A mobility edge in the imaginary axis in non-Hermitian Hamiltonians

Searching for Anderson localization of light in three dimensions has challenged experimental and theoretical research for the last decades. Here the problem is analyzed through large-scale numerical simulations, using a radiative Hamiltonian, i.e., a non-Hermitian long-range hopping Hamiltonian, well suited to model light-matter interaction in cold atomic clouds. Light interaction in atomic clouds is considered in the presence of positional and diagonal disorder. Due to the interplay of disorder and cooperative effects (sub- and super-radiance) a novel type of localization transition is shown to emerge, differing in several aspects from standard localization transitions which occur along the real energy axis. The localization transition discussed here is characterized by a mobility edge along the imaginary energy axis of the eigenvalues which is mostly independent of the real energy value of the eigenmodes. Differently from usual mobility edges it separates extended states from hybrid localized states and it manifests itself in the large moments of the participation ratio of the eigenstates. Our prediction of a mobility edge in the imaginary axis, i.e., depending on the eigenmode lifetime, paves the way to achieve control both in the time and space domains of open quantum systems.

Introduction.The interplay of opening and disorder in systems described by non-Hermitian Hamiltonians has been at the center of interest in many research fields, showing that non-Hermiticity can strongly affect the response of a system to disorder, inducing many counterintuitive effects [1][2][3][4][5][6][7][8][9][10][11].On the other side, Anderson localization [12] has been a beacon to understand closed disordered systems and has been at the focus of an ever increasing research community, ranging from condensed matter to acoustics, optics, and ultra-cold matter waves as well as quantum memories based on cold atoms [13][14][15][16][17][18][19][20][21].In the standard Anderson localization problem, an excitation can tunnel to nearest-neighbor sites placed in a regular lattice with disordered on-site energies (diagonal disorder).Depending on the value of the disorder strength, a mobility edge can be present at a specific energy: below this energy the eigenstates are localized, while above they are delocalized.
Extending the concepts developed for Anderson localization to open quantum systems still remains a challenge.Light has been an obvious candidate to study Anderson localization of non-interacting waves, which has triggered continuous efforts since the mid-80s [23][24][25][26][27][28][29][30][31][32][33].So far, Anderson localization of light in three dimensions however has resisted experimental observation.It has now been shown that pioneering experiments on Anderson localization of light [26][27][28] do not provide a signature for the Anderson transition in three dimensions [29][30][31][32][33] and the mere existence of an Anderson phase transition for light had even been questioned [34,35].Localization of light indeed presents many features which strongly differ from the standard Anderson localization of closed systems: i) in typical samples, scatterers have random positions in a three-dimensional volume, leading to positional disorder, ii) light induces complex long-range hopping between the sites, which in the case of two-level systems as scattering medium can lead to cooperative effects such as Dicke sub-and superradiance [36][37][38][39][40], iii) the excitation can escape from the system by photon emission, thus placing the problem of localization of light within the framework of open quantum systems.Both the longrange nature of the hopping and the opening can strongly affect the interplay of disorder and transport.Thus, the possibility to have a transition to localization in such systems is highly non-trivial.Specifically, cooperativity can affect the response of the system to disorder in a drastic way: while superradiant states show robustness to disorder [9,10], in the subradiant subspace long-range FIG. 1. (Color online) Localized Subradiant Eigenstates.A typical localized subradiant state for the scalar model and for Γ = 0.094Γ0, E = 0.1Γ0 and a participation ratio P R2 ≈ 7 for the case W/(Γ0b0) = 0.4 and N = 6400, ρλ 3 = 5 so that b0 ≈ 17.3.Upper panel: Three-dimensional representation of a localized eigenstate.The radius representing each atom is proportional to its excitation probability |Ψj(r)| 2 , also coded in color [22].Lower panel: A localized eigenstate projected on the x − y plane.
interaction is effectively shielded [41,42] and signatures of localization can emerge [6,7,41].In this letter, we shed new light on the problem of light localization in resonant scattering media by combining the positional disorder studied so far, with the initial ingredient of diagonal (on-site) disorder as considered by Anderson [12].With the aid of large scale simulations of up to 50000 atoms relying on a well-known radiative non-Hermitian Hamiltonian [43] able to take the vectorial nature of light into account, we show that the playground of the problem of light localization lies in the complex eigenvalues of the radiative Hamiltonian.Specifically, by analyzing the generalized participation ratio (GPR) of the eigenstates, a widely used figure of merit to study localization transitions, we show that a drastic transition in the behaviour of the GPR occurs at a specific imaginary energy value (decay widths) of the eigenstates.Such transition shares many analogies with what happens at real energy mobility edges and thus reveals the presence of a mobility edge along the imaginary axis in such systems.
The Model.We model light scattering in a 3D cold atomic cloud by considering N atoms randomly distributed inside a cube of volume V = L 3 , with a spatial density ρ = N/L 3 .When considering the interaction of atoms with the electromagnetic field, the full vectorial character of light should be taken into account.We will focus on atoms driven on a s → p dipole radiation transition which can be characterized by three degenerate levels in the excited state (labelled as α = x, y, z), each with a transition dipole moment (TDM) equal in coupling strength and perpendicular to the others [35].Thus we model each atom as a four-level system, with a ground state |g⟩ and three degenerate excited states |x⟩, |y⟩ and |z⟩.Corresponding TDM matrix elements are ⟨α| ⃗ µ|g⟩ = µê α , with α = x, y, z and the Cartesian unit vectors defined as êα .The radiative hamiltonian H vec describing dipole-dipole coupling in the single excitation approximation (see also [43]) is where E n,α are the atomic transition energies and Γ 0 is the radiative decay rate of a single atom.In Eq. ( 1), |n, α⟩ represents a quantum state where the nth atom is excited in its αth state, while all the other atoms are in their ground state.Interaction terms are non-Hermitian, namely In Eq. ( 2), k 0 = 2π/λ = E 0 /(hc) is the transition wavenumber (where λ is the wavelength of the atomic transition and E 0 being the average single atomic transition energy E 0 = ⟨E n,α ⟩, where the average is taken over disorder realization.r m,n is the distance between the mth and nth atom and rm,n is the unit vector joining them.Together with the vectorial model we also consider the scalar model [44].Even though the latter approximation neglects polarisation effects, it is appropriate in the dilute limit, where inter-atomic distances are larger than the optical wavelength λ, making near-field terms decaying as 1/r 3 negligible [22].The effective Hamiltonian which governs the interaction of the atoms with the electromagnetic field in the scalar approximation is characterized by complex long-range hopping terms V m,n decreasing as 1/r m,n with the distance, where the state |n⟩ stands for the n−atom in the excited state and all the other atoms being in the ground state, while is the interaction between the atoms at distance r m,n .The model in Eq. ( 3), known as the scalar model, has been introduced first by Foldy [45] and it has been used in several papers to describe cold atomic clouds in the dilute limit [46].Note that the vectorial model Hamiltonian has dimension 3N × 3N contrary to the scalar case which has dimension N ×N .Thus the scalar model allows us to investigate much larger system sizes with better statistics.For both models we can define the resonant mean free path l = 1/ρσ 0 (in the independent scattering approximation), where σ 0 = 4π/k 2 0 is the resonant scattering cross section in a simplified scalar model.Finally, we define the resonant optical thickness, b 0 , as the ratio between the system size L and the mean free path l.For the scalar model we have while for the vectorial model the resonant scattering cross section is σ 0 = 6π/k 2 0 and the optical thickness thus has to be corrected with respect to the scalar case : b Note that H contains both real and imaginary parts, which takes into account that the excitation is not conserved since it can leave the system by outgoing radiation.Its complex eigenvalues E = E − iΓ/2 describe the energies and line-widths (decay rates) of the eigenmodes of the system.We stress that even in the dilute limit ρλ 3 ≪ 1 we can have cooperative behaviour in the large sample limit (L ≫ λ), provided that the cooperativity parameter is b 0 ≫ 1.In this regime cooperative effects such as single-excitation sub-and superradiance become relevant [38,46,47].
In addition to the positional disorder of the atoms as studied previously [34,35], we now introduce an additional random diagonal disorder term in the Hamiltonian, which shifts the excitation energy of the atoms around its avarege value E 0 .Such diagonal disorder terms have not been given sufficient consideration in the context of localization of light, as engineering such effects is difficult in typical condensed-matter samples.However in cold atomic clouds, such on-site disorder can be realized by applying a speckle field coupling the excited state to an auxiliary excited state with convenient detuning, inducing thus random light shifts of the atomic resonances without inducing dipole forces in the ground state.Following the approach of the Anderson model on a lattice, we allow the site energies to fluctuate in the range of [−W/2, +W/2], where W is the strength of disorder.Ensemble averaging thus includes different realizations of the random position of the atoms and of site disorder.Within this model, we study both the eigenvalues as has been done in [34,35] as well as the eigenstates [48].
Localized Subradiant States.A striking illustration of the existence of localized states is given in Fig. 1, where we represent a typical localized subradiant eigenstate for the scalar model.The upper panel of Fig. 1 shows a 3D representation of a typical localized eigenstate, while the lower panel of Fig. 1 shows the projection of the squared wavefunction |ψ(r)| 2 on the x − y plane.While for zero diagonal disorder the vast majority of the states, which can be both superradiant or subradiant, are fully delocalized [22] for the spatial density considered, adding sufficient diagonal disorder leads to localization of the longerlived subradiant states.We observe that the localized peak, shown in the lower panel of Fig. 1, comes hand in hand with an extended tail, thus exhibiting a hybrid character, in agreement with Refs.[6,7,41].We note that the presence of such extended tails might strongly affect transport properties [49], for instance suppressing the exponential decay of transmission with the system size.Here we focus on the structure of the eigenmodes, leaving the analysis of the transport properties of subradiant localized states for a future work.
Mobility edge in the imaginary axis.In open systems, standard approaches to study localization such as the Thouless parameter should be applied with care [50].We therefore analyze the properties of the generalized participation ratio (GPR) of the eigenfunction ψ of the system [51,52], For localized eigenfunctions P R q is independent of the system size for all q, while in the delocalized regime P R q ∝ N q−1 .On the other hand, at the localization transition the GPR diverges with N as where d is the embedding dimension and D q defines the fractal dimension.Moreover, the distribution P (P R q /P R typ q ), where P R typ q = exp ⟨ln(P R q )⟩, is invariant at criticality in the large system size limit.This implies that the variance of the distribution P (ln(P R q )) is independent of the size at criticality [53][54][55][56][57][58][59][60], allowing for a precise identification of the critical point.
In order to have a general view of the localization properties of the eigenmodes of our system, we computed the GPR of all the eigenmodes for a specific disorder strength, and we plotted them as a function of their complex eigenvalues (real and imaginary part).A typical example of this analysis can be seen in Fig. 2, which shows a strong dependence of the P R 2 of the eigenmodes on the imaginary part of their eigenvalues, while the dependence on the real part is weak.Specifically we observe that the smaller is their imaginary part, the more the eigenmodes are localized.Note that the results of Fig. 2 refer to the scalar model, and a similar figure for the vectorial model can be found in [22].These results are consistent with previous findings about the interplay of super-and subradiance with disorder [6,7,9,10]: subradiant states are the ones most affected by disorder.The most interesting feature of this non-uniform response of the eigenmodes to disorder can be seen if one analyzes the typical value of P R q (P R typ q ) as a function of the decay widths.Since the optical thickness b 0 sets a relevant energy scale of the system (i.e. the spectral energy broadening prior to adding the diagonal disorder) [35], we considered different systems at a constant density and for a fixed value of the ratio W/(b 0 Γ 0 ).The results are shown in Fig. 3  the presence of a transition in the behaviour of the GPR: while the typical P R q of the eigenmodes is independent of the system size below Γ cr if W/b 0 is kept fixed (see vertical dashed line), it increases with the system size above Γ cr .These are precisely the same features present when analyzing the GPR of the 3D Anderson model (or other models displaying a localization transition) in correspondence of a mobility edge in the real energy.Thus our results points to the existence of a "mobility edge" in the imaginary axis.We checked that the imaginary mobility edge is independent of E around the band center, as shown in Fig. 2 and further discussed in [22].
We note that in the large density limit the results shown in Fig. 3(c) are extremely interesting, since they indicate that in presence of diagonal disorder, a localization transition can exist even in the large density limit for the vectorial case in absence of any magnetic field.This is at variance with what has been stated in [34] where no diagonal disorder was considered.Moreover the mobility edge in the imaginary axis, even in the large density limit, is well captured by Eq. (7).
In order to identify the critical decay width corresponding to the imaginary mobility edge, we performed a systematic analysis of the variance of the GPR vs the disorder strength W/Γ 0 for different densities, system sizes and ranges of decay widths.The variance of ln(P R q ) has been used in the literature to pinpoint the localization transition and it has been shown that, at the localization transition, the variance of ln(P R q ) is independent of the system size due to a universal distribution of the GPR [55].Similarly to Ref. [60], we use the crossing of rms(ln P R q ) close to its maximal value to locate the localization transition, see insets in Fig. 3.This allowed us to identify a critical decay width Γ cr .
We studied the critical decay widths as a function of disorder for different densities.The results are shown in Fig. 4(a).By fitting the numerical results we obtained an expression for the critical decay width: We note that the above expression cannot be extrapolated at small values of disorder since in that case the landscape of the GPR can only be understood analysing the whole complex plane [22].
We have also analyzed the GPR for different q values: q = 0.1, 0.6, 2, 5, 10 [22].For q ≥ 2 we always find a clear signature of a localization transition at a critical decay width, while for small values of q a localization transition is not observed.This reflects the hybrid nature of the localized eigenstates: indeed together with a localized peak, an extended tail is present.The GPR for large values of q is more sensitive to large values of |ψ| 2 , thus it describes the behaviour of the localized peak, whereas the GPR for small values of q is sensitive to small values of the wave function amplitudes and thus to the wave function tails.Since the tails are always extended (delocalized), no localization transition is seen for small q [22].In order to further confirm the above picture, we have computed the fractal dimension D q as a function of the decay widths.The results are shown in Fig. 4(b).As one can see for q ≥ 2 a transition in the fractal dimension is seen from zero to a value larger than one, while for q < 1 no transition is observed, confirming the hybrid nature of the eigenmodes of the system.Note that in the extended phase, even for q = 2, 5, D q is different from d = 3 indicating that the wave function are never fully extended.In other words, the eigenfunctions are always multifractal both below and above criticality: Γ cr marks the transition from a frozen phase (where the the GP R is independent of N for sufficiently large q), to a weakly multifractal phase (with a narrow distribution of fractal dimensions D q ) [61].
Conclusions.We considered a well known radiative non-Hermitian Hamiltonian model to describe coherent multiple scattering of light in cold atomic clouds at low excitation level.Our results give new insights on the problem of localization in open quantum systems under the interplay of non-Hermiticity and disorder.A novel kind of localization transition has been identified, occurring at a critical lifetime (or inverse decay rate Γ) of the eigenmodes of the system, i.e. along the imaginary energy axis.A single-parameter scaling was found for the critical decay rate Γ cr /Γ 0 [Eq.(7)] for the localization transition, which is given by W/(b 0 Γ 0 ), in contradiction to what could be expected from ρ or b 0 separately.The localization transition identified here in a realistic model of light matter interaction shares many analogies with the Anderson transition in 3D lattices and with localization transitions in long-range interacting systems, such as in the power-banded random matrix model [55,61], but also important differences: the localization transition is signalled by the behaviour of the GPR for large q values (larger or equal than 2) and not for small q values (less than 2).We attribute this feature to the fact that the eigenmodes are not fully localized but that have a hybrid character, with a localized peak and an extended tail.A precise characterization of the shape of these eigenmodes will be the topic of a future work.Despite these differences, our results indicate the existence of a novel kind of localization transition occurring along the imaginary energy axis which is independent of the real energy (around the band center) for sufficiently large values of diagonal disorder and optical thickness.The existence of a mobility edge in the imaginary axis found in this Letter certainly constitutes a novel feature in the field of localization in open quantum systems.Further research will be necessary to assess the impact of our results.For instance the general conditions for this mobility edge in the imaginary axis to arise in open quantum systems should be investigated both in the single excitation and many excitation regime and for different topology and dimensions and critical exponents determined.Our findings are relevant not only from a fundamental point of view but also for applications, e.g. to achieve efficient energy storage, quantum memory, quantum simulation and sensing devices.

I. EXTENDED SUBRADIANT STATE
Here we show an example of a typical extended subradiant state for the scalar model in absence of diagonal disorder [W/(b 0 Γ 0 ) = 0], see Fig. S1.This figure should be compared with Fig. 1 of the main text where a typical localized subradiant state with W/(b 0 Γ 0 ) = 0.4 is shown.Comparing the two figures one can see that disorder in the transition frequencies of the atoms can induce localized states in the subradiant subspace.Both in Fig. S1 of this supplementary material and Fig. 1 of the main text, in the upper panels each atom is shown by a small sphere.The probability |Ψ j (r)| 2 for the eigenstate to be on that atom is given by the color and the radius R of the sphere according to the relation , where |Ψ j (r)| 2 max is the maximal probability for the case W/(b 0 Γ 0 ) = 0.4.This normalization relation was chosen to improve visibility.In the lower panels the projection on the x − y plane of |Ψ j (r)| 2 on a grid of 60×60 is shown.To improve the quality of the representation, each grid point has been averaged by the surrounding points, with a weighting inversely proportional to their distances squared.

II. MOBILITY EDGE IN THE IMAGINARY AXIS: SCALAR MODEL
In order to analyze the localization transition we consider the typical value of the Generalized Participation Ratio (GPR) of the eigenmodes of the system, see Eq. ( 5) in the main text.We note that the eigenfunctions of the non-Hermitian Hamiltonian represent the projection of the total eigenfunctions on the single excitation manifold of the atomic degrees of freedom.Thus the quantity |ψ k | 2 which is used to compute the P R q represents the conditional probability to find the system on atom k, given that one quantum of excitation is stored in the system.The state |k⟩ is the state where the atom k is excited while all the other atoms are in the ground state.
In order to have a general view of the localization properties of the eigenmodes of our system, we computed the GPR of all the eigenmodes for a specific value of the disorder, and we plotted them as a function of their complex eigenvalues (real and imaginary part).A typical example of this analysis can be seen in Fig. S2(a), which shows a strong dependence of the P R 2 of the eigenmodes on the imaginary part of their eigenvalues, while the dependence on the real part is weak.Specifically we observe that the smaller is their imaginary part, the more the eigenmodes are localized.We also analyzed the typical value of P R q (P R typ q ) as a function of the decay widths for different systems at a constant density and for a fixed value of the ratio W/(b 0 Γ 0 ).The results are shown in Fig. S2(b,c): they clearly indicate a localization transition.While the typical P R q of the eigenmodes is independent of the system size below Γ cr if W/(b 0 Γ 0 ) is kept fixed (see vertical dashed line), it increases with the system size above Γ cr .In order to identify the critical decay width corresponding to the imaginary mobility edge, we performed a systematic analysis of the variance of the GPR vs. the disorder strength W/Γ 0 for different densities, system sizes and ranges of decay widths.The variance of ln(P R q ) can be used, see main text, to pinpoint the localization transition: we use the crossing of rms(ln P R q ) close to its maximal value to locate the localization transition, see insets in Fig. S2(d).This allowed us to identify a critical decay width Γ cr , see details in the main text.Using Eq. ( 5) in the main text and performing a scaling analysis we can determine the fractal dimension as a function of the decay widths for different q values: D q = d q−1 ln(P R q )/ ln(N ).
The rescaled typical GPR P R q /N Dq(q−1)/3 , where D q is the fractal dimension computed at criticality, is shown in the insets of Fig. S2(b,c).As one can see the re-scaled P R q nicely cross at the critical decay width.
In Fig. S3 the typical value of P R q is shown for different values of q for the case ρλ 3 = 5, W/(b 0 Γ 0 ) = 0.3.As one can see a clear signature of a localized-delocalized transition is shown for large values of q = 2, 5 (lower two panels), while for small values of q = 0.1, 0.6, P R q increases with the system size for all decays widths.As discussed in the main text, in the open quantum system considered here, localized eigenmodes have a hybrid nature, with a localized peak and an extended tail.Small q values are sensitive to the tails and thus reveal the extended character of the eigenmodes, while large q values are more sensitive to the peak, thus revealing the localized character of the eigenmodes.Note that even if the P R q is never independent of the system size for small q values, their dependence on the decay widths has a change of slope in correspondence of the imaginary mobility edge, see vertical dashed lines in the upper panels in Fig. S3.The mobility edge in the imaginary axis can be also analyzed considering the generalized inverse participation ratio (GIPR) of the eigenfunction ψ of the system, For localized eigenfunctions IP R q is independent of the system size for all q, while in the delocalized regime IP R q ∝ N 1−q .In Fig. S4(a,b,c) the typical IPR is shown for q = 2, 5, 10, showing the mobility edge in the imaginary axis with the same critical decay width as computed in Eq. ( 7) in the main text, see vertical dashed line.In Fig. S4(d) the root mean square of ln(IP R q ) is shown for different q values.As one can see our estimation for the critical decay width (see vertical dashed line) indicates fairly well the crossing of rms(ln(IP R q )) for different N even if a slight difference between the crossing for q = 2 and q = 5, 10 is visible.A detailed investigation of this effect is beyond the scope of this manuscript and it will be investigated in future work.Note that the root mean square of ln(IP R q ) and ln(P R q ) are the same since the inverse participation ratio and the participation ratio are just the inverse of each other.
Finally, it is important to note that our estimation of the critical decay width as a function of the disorder strength [Eq.(7) in the main text] cannot be extrapolated to small values of disorder.Indeed for very small disorder the mobility edge in the imaginary axis is not defined, see Fig. S5.In general the playground for the localization of open quantum systems is the complex plane, see Fig. S5, where the typical P R q=2 for the case of zero disorder is shown in the complex plane.As one can see comparing this figure with Fig. 2 of the main text and with Fig. S2(a), for zero diagonal disorder no clear mobility edge in the imaginary axis is present.On the other side, one cannot exclude the presence of other mobility edges along different boundaries in the complex plane (this topic is outside the focus of the current manuscript and it will be investigated in a future work).

III. MOBILITY EDGE IN THE IMAGINARY AXIS: VECTORIAL MODEL
When considering the interaction of atoms with the electromagnetic field, the scalar model, see Eq. (3) in the main text, is valid for dilute systems.In general the full vectorial character of light should be taken into account.
Here we describe in detail the radiative hamiltonian H vec describing light-matter interaction in atomic-like systems for weak fluence (single excitation approximation), see also [1].We will focus on atoms driven on a s → p dipole radiation transition which can be characterized by three degenerate levels in the excited state (labelled as α = x, y, z), each with a transition dipole moment (TDM) equal in coupling strength and perpendicular to the others [2].
Thus we model each atom as a four-level system, with a ground state |g⟩ and three degenerate excited states |x⟩, |y⟩ and |z⟩.Corresponding TDM matrix elements are α ⃗ µ g = µê α , with α = x, y, z and the Cartesian unit vectors defined as êα .
To model cold atomic clouds, here we consider an ensemble of atoms randomly placed in a 3D box (positional disorder).Note that for the vectorial model the optical thickness has to be modified with respect to the scalar case, and we have: b , where the optical thickness for the scalar case is given in Eq. ( 4) of the main text.The Hamiltonian which takes the vectorial nature of light into account, describing an ensemble of atoms interacting with light can be written as [2]: In Eq. (S3), k 0 = E 0 /(hc) is the transition wavenumber (with E 0 being the mean atomic transition energy E 0 = ⟨E n,α ⟩), r m,n is the distance between the mth and nth atom and rm,n is the unit vector joining them.Note that this Hamiltonian has dimension 3N × 3N contrary to the scalar case which had dimension N × N .Using Eq. (S2) we have analyzed an ensemble of atoms in a 3D box occupying random positions.Diagonal disorder has been also considered allowing the energies E n,α to fluctuate in the range of [−W/2, +W/2], where W is the strength of disorder.The typical values of the generalized participation ratio has been computed by computing the probability of the excitation to be on every atom.Fig. S6 shows the results for the vectorial case ρλ 3 = 5, W/(b 0 Γ 0 ) = 0.8.As one can see, clear signatures of the mobility edge in the imaginary axis are shown for the typical value of P R q=2,5 and for rms(ln P R q=2,5 ).Note that our estimation of the critical decay width corresponding to the imaginary mobility edge given in Eq. (7) in the main text is in very good agreement with the numerical data, see vertical dashed lines in Fig. S6.Now we turn our attention to the large density case, where the full vectorial model is particularly relevant since the scalar model is a good approximation only in the dilute limit (small densities).It has been claimed that, in absence of a magnetic field, localization is not possible in the vectorial model of light [3].Nevertheless, here we show that, even in the large density limit, the introduction of diagonal disorder induces a mobility edge in the imaginary axis which is well captured by Eq. ( 7) in the main text, see Figs.S7,S8.
On the other hand, in absence of disorder and for large densities, the localized features of the eigenmodes of the system are more complicated to capture.Again the real playground is the complex plane, see Fig. S9.The possibility to find localized states even in absence of diagonal disorder, and in presence of positional disorder only, will be investigated elsewhere.

IV. ON THE NATURE OF THE MOBILITY EDGE IN THE IMAGINARY AXIS.
Usually in open Anderson models [4], the excitation can escape the system only from the boundaries, so that the decay widths are proportional to the probability of a state to be on the boundaries.As a consequence of this, most of the localized states also have very long lifetimes  Note that E k is the difference between the real part of the eigenvalues and E0.
(similar to subradiant states), since their probability to be on the boundaries is exponentially small.On the other side, the model studied here, see also Ref. [5], strongly differs from the previously studied models of localization in open systems, since in our case the excitation can escape from any site and not only from the boundaries.For instance in our model a fully localized state on one site has a decay width equal to Γ 0 , independent of the system  Note that E k is the difference between the real part of the eigenvalues and E0.

size.
In order to clarify the difference between our model and previously studied open 3D Anderson models, let us consider a 3D cubic Anderson model with leads connected to one of its side as in Ref. [4].Let us assume that the disorder is such to create a mobility edge at energy E c .Clearly the decay width of the states will be very small for E < E c , while they will be large for energy E > E c , and correspondingly a mobility edge could also be found in the imaginary axis if one plots the participation ratio P R 2 vs the decay widths.But in this case to use E or Γ is just a different way to label the states.On the other side our mobility edge has a completely different nature since it is independent of the real energy of the states in a wide energy range around the energy center, and it only depends on their imaginary energy, i.e. the lifetime of the eigenmodes of the open system.
We note that the dependence of the P R 2 on the lifetime of the subradiant eigenmodes is a novel feature, which has not been captured by the toy model of Refs.[6,7].Indeed, in the open 1D and 3D Anderson model analyzed in Refs.[6,7], the sub-and superradiant modes were segregated in two regions, whereas in the present case, no gap between sub-and superradiant modes exists.
We also note that in the closed Anderson 3D model, the P R 2 diverges at a finite energy corresponding to the mobility edge.In our case, the P R 2 diverges at a finite decay width (corresponding to the imaginary part of the complex eigenvalues of the system), thus we use the term mobility edge in the imaginary axis in analogy with the localization transition in closed systems which occurs along the real axis.In the case of a closed system, such as the standard Anderson model, the behavior of the P R 2 reflects the transport properties of a system in a direct way: when the P R 2 increases with the system size, transmission will be diffusive or ballistic, while if the P R 2 is independent of the system size, transmission is exponentially suppressed with the system size due to localization.In the case of open systems, described by a non-Hermitian Hamiltonian, the P R 2 has a more indirect link to transport properties since the eigenmodes are not fully localized but they have an hybrid nature as discussed in the main text.A discussion of the transport properties of hybrid states can be found in Ref. [9].We do not aim to discuss this further in this manuscript, we just note that 1/(E − H) is the propagator for the excitation in the system.For this reason, a change in the structure of the eigenmodes of H as signaled by the P R 2 represents a real physical change in the way excitations propagate through the system.
For the model considered here, the increase of the P R 2 with Γ might be explained by the increase of the mean level spacing of the eigenvalues with Γ, see Fig. S10.Note that in Fig. S10 we compute the mean level spacing in the complex plane of the complex eigenvalues of the system.Indeed perturbation theory in the case of non-Hermitian Hamiltonian shows that it is the distance in the complex plane which determines the strength of perturbations [8,10].Moreover we checked that also the mean level spacing in the real axis increases with the decay width.The increase of the mean level spacing with the decay width is due to the fact that superradiant states have a stronger coupling to the photon field, so that their energy spreads much more than subradiant states, which are partially shielded from the interaction [11].Thus for a fixed amount of disorder, the states with lower Γ are more easily mixed by disorder than the states with a larger Γ.Nevertheless this argument cannot explain the emergence of a mobility edge in the imaginary axis.More work is needed to deepen the understanding of this novel feature and to understand the general requirements for the mobility edge in the imaginary axis to emerge.
We also note that a highy non-uniform mean level spacing is typical for systems with long-range interactions.For instance even a finite energy gap can be induced in such systems [10,11].Localization even in presence of long-range interactions has been discussed for subradiant states in [6,7] and in general, in the framework of a shielding effect in [11,12] and more recently in [13].
Finally we would like to point out that often in literature, localization properties are studied by means of the Thouless parameter [3,14].The Thouless parameter requires only the eigevalues of the system, nevertheless it should be used with care in open systems.Indeed, in presence of absorption or other sources of leakage, the decay width of the states should be properly redefined in order to take into account only the leakage from the boundaries, see discussion in Ref. [15].Thus, analyzing the structure of the eigenmodes, as we did here, is a much reliable mean to study localization in open quantum systems.The mean level spacing in the complex plane D is plotted vs the decay widths Γ (red circle) for the case N = 3200, ρλ 3 = 5 and W = 0.The mean level spacing has been computed by counting the number of complex eigenvalues per unit area in the complex plane for −0.1 < (E −E0)/Γ0 < 0.25 and different ranges of Γ.The mean level spacing D has been obtained by taking the square root of the inverse density of complex eigenvalues.Note that an increase of the mean level spacing is observed even if one computes the distance in the real energy axis of the complex eigenvalues.

FIG. 2 .
FIG. 2. (Color online)Participation ratio in the complex plane: The participation ratio P R = P Rq=2 of the eigenstates for the scalar model is shown in the complex plane (E k , Γ) of the complex eigenvalues for N = 283 , ρλ 3 = 5, b0 ≈ 26.06 and W/(b0Γ0) = 0.3.Note that the horizontal energy axis is shifted by E0.The critical width for the transition to localization [Eq.(7)] is indicated by the red horizontal line.Note that E k is the difference between the real part of the eigenvalues and E0.

FIG. 3 .
FIG. 3.(Color online) Mobility edge in the imaginary axis.Typical GPR for q = 2 (panel a,b,c) and q = 5, 10 (panel d) as a function of the decay width of the eigenstates are shown both for the vectorial (a,c) and scalar model (b,d) for different number N of atoms and constant density, see legend.The vertical black dashed line indicates the critical width [Eq.(7)].In the insets the root mean square of ln(P Rq) is shown as a function of the decay width of the eigenstates.In panels (a,b,d), for each N the eigenvalues in the region −b0/4 < (E − E0)/Γ0 < b0/4 were considered, while in panel c)the region −7 − b0/4 < (E − E0)/Γ0 < 7 + b0/4 were considered.
, both for the vectorial and the scalar model, for different system sizes at constant density.The results clearly indicate

FIG. 4 .
FIG. 4. (Color online) Panel a), critical decay width: the critical decay width for the localization transition is shown as a function of the normalized disorder strength for different densities and for both the vectorial and the scalar model (see legend).The precision with which we determined the critical decay width is always below ±0.015.Panel b): fractal dimension: the fractal dimension as a function of the normalized decay width is shown for the case ρλ 3 = 5, W/(b0Γ0) = 0.5.The fractal dimension has been extracted from the size dependence of the GPR for different values of q.The vertical black dashed line indicates the critical width [Eq.(7)].

2 P
FIG. S1. (Color online) Representations of a typical extended subradiant state (for the scalar model) similar to Fig. 1 in the main text.Upper panel: Three-dimensional representation of an eigenstate.The radius representing each atom is proportional to its excitation probability |Ψj(r)| 2 .Lower panel: An eigenstate projected on the x−y plane.Here N = 6400, ρλ 3 = 5 so that b0 ≈ 17.3 and W/(b0Γ0) = 0.For the state shown in both panels we have E/Γ0 = −0.0758,Γ/Γ0 = 0.05.The participation ratio P R2, defined in Eq. (5) of the main text, of the state shown in this figure is P R2 = 1941.

3 FIG
FIG. S2. (Color online)Mobility edge in imaginary axis.(a) Participation ratio P R = P Rq=2 of the eigenstates in the complex plane (E, Γ) of the eigenvalues of each state for N = 283 , ρλ 3 = 5, b0 ≈ 26.06 and W/(b0Γ0) = 0.5.Note that E k is the difference between the real part of the eigenvalues and E0.The critical width for the transition to localization [Eq.(7) in the main text] is indicated by the red horizontal line.(b,c) Typical GPR for q = 2 (panel b) and q = 5 (panel c) as a function of the decay width of the eigenstates for W/(b0Γ0) = 0.5, ρλ 3 = 5.The vertical black dashed line indicates the critical width [Eq.(7)  in the main text].In the insets the P Rq/N Dq (q−1)/d , where Dq is the fractal dimension computed at criticality, is shown in the region around the localized-delocalized transition.(d) Participation ratio fluctuations.The root mean square of ln(P R2) is shown as a function of the decay width of the eigenstates.The inset shows an enlargement of the same panel around the transition.All panels refer to the case W/(b0Γ0) = 0.5, ρλ 3 = 5 and different number N of atoms, see legend.In panels (b,c,d), for each N the eigenvalues in the region −b0/8 < (E − E0)/Γ0 < b0/8 were considered.

FIG
FIG. S5. (Color online) Absence of a mobility edge in the imaginary axis in absence of disorder for the scalar case.Participation ratio P R = P Rq=2 of the eigenstates (see legend on the right) in the complex plane (E/Γ0, Γ/Γ0) of the eigenvalues of each state for N = 32 3 , ρλ 3 ≈ 5.05, b0 ≈ 30 and W = 0. Note that E k is the difference between the real part of the eigenvalues and E0.
FIG. S6. (Color online) Mobility edge in the imaginary axis for the vectorial case.Panel (a): The root mean square of ln P Rq for q = 2, 5 is shown as a function of the decay width of the eigenstates for W/(b0Γ0) = 0.8, ρλ 3 = 5.Panel b,c): Generalized Participation ratio for q = 2 (panel b) and q = 5 (panel c) as a function of the decay width of the eigenstates for W/(b0Γ0) = 0.8, ρλ 3 = 5.Here the typical P Rq is averaged over the range −b0/4 < (E − E0)/Γ0 < b0/4.The vertical black dashed line in all panels indicates the critical width obtained from Eq. (7) in the main text.Data in all panels refer to the vectorial model, see Eq. (S2).
FIG. S7. (Color online) Mobility edge in the imaginary axis for the vectorial case.Participation ratio P R = P R2 of the eigenstates (see legend on the right) in the complex plane (E/Γ0, Γ/Γ0) of the eigenvalues of each state for N = 25 3 , ρλ 3 = 40, and W/(b0Γ0) = 0.5.The horizontal red line indicated the critical decay width, see Eq. (7) in the main text.Note that E k is the difference between the real part of the eigenvalues and E0.
FIG. S8. (Color online) Mobility edge in the imaginary axis for the vectorial case.Panel (a,b): Typical generalized participation ratio for q = 2 (panel a) and q = 5 (panel b) as a function of the decay width of the eigenstates.Panel (c,d): The root mean square of ln P Rq for q = 2, 5 is shown as a function of the decay width of the eigenstates.Here the typical P Rq is averaged over the range −7 − b0/4 < (E − E0)/Γ0 < 7 + b0/4.All data refer to the case W/(b0Γ0) = 0.5, ρλ 3 = 40.The vertical black dashed line in all panels indicates the critical width obtained from Eq. (7) in the main text.Data in all panels refer to the vectorial model, see Eq. (S2).
FIG. S9. (Color online) Absence of a mobility edge in the imaginary axis in absence of disorder for the vectorial case.Participation ratio P R = P Rq=2 of the eigenstates (see legend on the right) in the complex plane E/Γ0, Γ/Γ0 of the eigenvalues of each state for N = 25 3 , ρλ 3 = 40, and W = 0.Note that E k is the difference between the real part of the eigenvalues and E0.

3 FIG
FIG. S10.(Color online)Increase of the mean level spacing with the widths.The mean level spacing in the complex plane D is plotted vs the decay widths Γ (red circle) for the case N = 3200, ρλ 3 = 5 and W = 0.The mean level spacing has been computed by counting the number of complex eigenvalues per unit area in the complex plane for −0.1 < (E −E0)/Γ0 < 0.25 and different ranges of Γ.The mean level spacing D has been obtained by taking the square root of the inverse density of complex eigenvalues.Note that an increase of the mean level spacing is observed even if one computes the distance in the real energy axis of the complex eigenvalues.
m̸ =n) α,β∈{x,y,z} V m,n,α,β |m, α⟩ ⟨n, β| , (S2) where E n,α are the atomic transition energies and Γ 0 is the radiative decay rate of a single atom.In Eq. (S2), |n, α⟩ represents a quantum state where the nth atom is excited in its αth state, while all the other atoms are in their ground state.Interaction terms are non-Hermitian, namely