Renormalization to localization without a small parameter

We study the wave function localization properties in a d-dimensional model of randomly spaced particles with isotropic hopping potential depending solely on Euclidean interparticle distances. Due to generality of this model usually called Euclidean random matrix model, it arises naturally in various physical contexts such as studies of vibrational modes, artificial atomic systems, liquids and glasses, ultracold gases and photon localization phenomena. We generalize the known Burin-Levitov renormalization group approach, formulate universal conditions sufficient for localization in such models and inspect a striking equivalence of the wave function spatial decay between Euclidean random matrices and translation-invariant long-range lattice models with a diagonal disorder. Copyright A. G. Kutlin and I. M. Khaymovich. This work is licensed under the Creative Commons Attribution 4.0 International License. Published by the SciPost Foundation. Received 23-01-2020 Accepted 25-03-2020 Published 01-04-2020 Check for updates doi:10.21468/SciPostPhys.8.4.049


Introduction
After more than six decades of successful and intense study of Anderson localization (AL) passed from Anderson's seminal work [1] this field still embodies many puzzles and unexpected surprises such as the correspondence between many-body localized (MBL) systems and AL models on hierarchical structures like a random regular graph (RRG), including the presence in both counterparts of a whole phase of the subdiffusive wavepacket spreading in the finite range of parameters [2][3][4][5][6], which is absent in single-particle models on finite-dimensional lattices, and hot debates about the presence of a putative extended phase violating ergodicity [7][8][9][10][11][12][13][14][15][16] with non-trivial multifractal wavefunctions, claimed in other papers to be just a finite-size effect [17][18][19][20][21][22][23][24][25][26][27] due to a critical regime [13] close to the localization transition. Another surprise 1 of recent years in AL community is the presence of robust localization of all bulk eigenstates in long-ranged (e.g., dipolar) systems [32] beyond the convergence of the standard locator expansion and of the resonance counting [1,33,34]. For lattice models with diagonal disorder this phenomenon is directly linked in the literature to the effects of cooperative shielding [35] 2 and the emergence of localization due to correlations in hopping [36,37] of these long-range models. However, this question for the models with off-diagonal disorder (e.g., due to disordered positions of lattice sites) is still open.
In this work, we address exactly this question via the discussion of the localization properties of systems described by Euclidean random matrices (ERM) [38], i.e. the systems of particles randomly distributed in d-dimensional space with the hopping of single-particle excitations, which solely depends on the Euclidean interparticle distance. Due to quite general description, ERMs cover a considerable class of physical models and arise naturally in various systems, e.g., in the ones with non-crystalline structures like gases, liquids, amorphous materials, and glasses. Although such models sometimes arise in the systems with shortrange interactions, such as elastic networks [39], jammed soft spheres [40] or magnetic vortex plasma [41], more commonly ERMs are used to describe the long-range models. Indeed, longrange ERMs are applied to the analysis of the systems of particles with Coulomb interactions in two-dimensional irregular confinement [42], disordered classical Heisenberg magnets with uniform antiferromagnetic interactions [43], systems with dipole-dipole interactions such as dipolar gases [44], systems of ultracold Rydberg atoms [45] and so on. Even the effects of photon localization in atomic gases [46] are described by a long-range ERM. Although the ERM model itself was introduced [38] back in 1999 and its spectral properties are studied quite deeply [47][48][49], the analysis of the wave function properties, including their spatial structure and localization, crucial for above mentioned applications to physical models is barely carried out and represented in the literature only by a couple of numerical works [32,50] or in quite restricted particular cases [51]. The present paper is aimed to fill this gap providing a generic analytical approach.
The problem with the analysis of the eigenstate structure in ERMs is caused by the absence of a small parameter. Indeed, unlike the models with the diagonal disorder, there is no way to treat the ERMs without the diagonal potential with the locator expansion approximation even for the infinitesimally small hopping term, due to the ideal resonance of all bare diagonal levels. This fact can be understood on the example of low-dimensional models with translation-invariant polynomial hopping which show localization of all bulk wavefunctions for any finite disorder either in the diagonal potential [32,35,36,52] or in the position of the lattice sites [32], whereas in the complete absence of the disorder these models are translation- 1 There are many other surprises such as emergence of multifractality in long-range static [28][29][30] or short-range driven [31] models with quasiperiodic potentials but we focus on the one relevant for our consideration. 2 In this case a top energy level keeps delocalized even at strong disorder due to its energy diverging with the system size and shields the rest levels from the hopping terms.
invariant and, hence, delocalized. To overcome the above mentioned principal difficulty we generalized the known renormalization group (RG) approach (developed by Levitov [33,34] and extended by Burin and Maksimov (BM) in [53] and by Mirlin and Evers in [54]) to the case of absence of a small parameter and to generic smooth Euclidean hopping term. Consequently, we show that for all ERMs with quite smooth potential 3 the bulk spectral eigenstates show localization.
The main idea behind this is similar to the one developed in [36], where the presence of (the measure zero of) delocalized states with energies diverging with the system size (not only the top energy level) at the either spectral edge gives the main contribution to the hopping term and is shown not to bring the system to the delocalization. In that case the effective hopping for the bulk spectral states can be obtained by the matrix-inversion trick developed in [36] which rewrites the eigenproblem in a special form, non-linear in eigenvalues, inverting all the high-energy contributions to the hopping term.
In the case of the current work on ERMs, the absence of a small parameter does not allow us to use the same technique and we have to develop a renormalization group (RG) approach. The resulting renormalization of the hopping terms is shown to evolve in such a way that the most of their spectral weight goes to the spectral edge states with energies increasing with the "renormalization scale" (system size) as in the matrix-inversion trick or the cooperative shielding. Thus, this significantly reduces the spectral weight of the hopping term in the bulk of the spectrum and localizes bulk spectral states.

Main idea
The cornerstone idea of the renormalization group approach with respect to AL in random matrix problems [33,34,53,54] is to rewrite the Hamiltonian of the system in such a form which is invariant under the iterative diagonalization procedure. The latter diagonalization procedure represents an elementary step of the renormalization scheme and thus should be done analytically as precise as possible. This often implies an exact diagonalization of certain 2 × 2 matrix blocks, which take into account most resonant levels hybridizing with the current one. This diagonalization procedure is crucially based on the assumption of isolated single resonant pairs of levels, which has a certain range of validity. The approximation can be formulated as follows: typically, for each iteration i and any energy level represented by a diagonal matrix element i n there is the only resonant level i m , m = n, such that the absolute value of the off-diagonal hopping element t i nm between them is comparable or larger than the interlevel spacing | i n − i m |. The RG procedure diagonalizing the initial problem is formed by a set of consecutive elementary diagonalizations of the resonant level pairs.
Further we consider a random matrix Hamiltonian of a general form with the deterministic real-valued function f (r) which depends only on the Euclidean distance r i j = |r i j | = |r i − r j | between some sites in d dimensions, indices i and j numerate all N ddimensional sites r i . The randomness in the model is given both by the off-diagonal elements through the positions of sites r i uniformly distributed in d-dimensional cube with the mean density equal to unity, and by the random bare on-site energies i with zero mean, dependent or independent of f (r) and r i . In particular, i could be even all equal to zero, as in the power-law Euclidean (PLE) model considered in [32] with d = 1 and f (r) = r −a . Unlike the models with translation-invariant hopping and only diagonal disorder [35][36][37]53], the above model does not necessarily have a small parameter, and, hence, the approximation of the single resonances does not necessarily applicable from the first steps of RG. In terms of the wave functions it means that the localized ones (if any) can have extended "heads" of finite size R 0 which will be determined later, where eigenstate do not decay, but may, e.g., oscillate. These heads have a complex internal structure which cannot be obtained within RG approach, because at such distances the approximation of single resonances may fail. To overcome this difficulty, we introduce the following preliminary step before employing RG: we rewrite the Hamiltonian in a form H = H 0 + V in such a way that, being expressed in the eigenbasis | 0 n 〉 of H 0 , it is invariant under the iterative diagonalization of resonant blocks and the approximation of single resonances is satisfied. This allows further RG treatment in the form of [33,53,54] and, thus, show the localization of the bulk of the states written in the basis | 0 n 〉. If, in addition, eigenvectors of H 0 are exponentially localized in the initial basis |c i 〉, one can equally consider the localization and the eigenstate spatial structure in either of bases | 0 n 〉 or |c i 〉. In this case one can forget about initial Hamiltonian and use the effective one instead.
To obtain the effective renormalizable Hamiltonian which satisfies all above mentioned conditions, we consider H 0 as the initial H cut at r i j ≤ R 0 and rewrite the original Hamiltonian in a form where R 0 is a cutoff radius at this zeroth step, and Since H 0 we used to obtain this form has short-range hopping in the original basis, the states | 0 n 〉 are assumed to be localized with the localization scale of the order of R 0 . Here we should note that even the worst case of the initial bare energies being strictly zero i = 0, corresponding to all bare sites being in perfect resonance already at the first RG step, is covered by this method. Indeed, taking R 0 = 1 one can easily diagonalize H 0 and get (i) exponentially localized eigenstates | 0 n 〉 and (ii) non-singular density of states (DOS) formed by nearly uncorrelated eigenvalues 0 n . Further we restrict our consideration to the most relevant case of smoothly varying hopping potentials f (r) at the scale R 0 4 and neglect the difference between r nm and r i j due to localized nature of the wavefunctions 〈 0 n |c i 〉 and 〈c j | 0 m 〉 at the cutoff radius R 0 r i j , r nm , see Fig. 1.
, and the effective Hamiltonian, Eq. (3), takes the form Here q 0 n = i ψ 0 n (i) is an effective "charge" of the state | 0 n 〉, with ψ 0 n (i) = 〈c i | 0 n 〉. As we show below, this Hamiltonian is renormalizable, with the effective charges being the renormalization parameters. The zeroth-step cutoff radius R 0 should be determined in such a way that the approximation of the single resonances is valid for the first step of the renormalization group.
First, we proceed to the renormalization group scheme which, from this point, is quite straightforward and leave the problem of the zeroth-step cutoff radius determination and the range of validity of the effective charge approximation for a further discussion. Assuming that on the ith iteration the renormalization group Hamiltonian H i has a form with R i+1 R i and the approximation of single resonances is valid, see diagonalization of the corresponding 2 × 2 resonant block coupling i n and i m levels This forms an elementary step of RG procedure which gives both the spectrum of the effective Hamiltonian H eff , Eq. (5), and the asymptotic form of tails of its localized eigenstates 5 . Indeed, for R i R 0 , when all the wavefunction heads are eventually formed, the strong resonances are rare and typical values of θ in (7c) are small. As a consequence, the typical wave functions transform as , the tails are determined by the effective hopping where is the squared effective charge for the state with energy and averaging (denoted by 〈. . .〉) is taken over disorder realizations and index n, and ν R i ( ) = 〈δ( − i n )〉 is the density of states (DOS) at ith RG step. Note that the energy dependence of the effective hopping is not accidental as there are few delocalized states at the spectral edge for which the RG approach is not applicable.
In determination of the spatial decay we should take into account the difference in wavefunction averaging. For the typical averaging ψ 2 n,typ (r m ) = exp 〈ln ψ 2 n (r m )〉 the eigenstate decays proportionally to t i eff ( ), (8) while for the mean averaging one have to take into account strong resonant contributions and obtain (due to Breit-Wigner-like profile of wavefunctions)

Basic equations
To determine the evolution of the effective charges we first write the equation for the probability P(q, ; R)dqd of a state at a certain RG step with the cutoff radius R to have energy and charge in the intervals ( , + d ) and (q, q + dq), respectively. Due to the hybridization (7) of resonant pairs the evolution of the probability distribution at one RG step takes the form Here, for brevity, we omit upper indices i and i + 1 and, instead of R i and R i+1 , write R 1 and R 2 . The integration by d d r nm is carried out over the whole region of the d-dimensional space in the interval R 1 < r nm < R 2 . Equation (12) provides an exact recipe to calculate the distribution function P(q, ; R) at all steps of the renormalization scheme provided the approximations of effective charges and single resonances are valid. Needless to say that due to this exactness and an overall complex structure of the equation, its analytical solution is extremely tough to obtain without further approximations. However, since the quantities of primary interest are few first moments of the distribution function and not the distribution function itself, we can, using Eq. (12), write similar equations for the moments and then try to solve them, exactly or approximately. For example, it can be directly seen from Eq. (12) that the average eigenenergy 〈 〉 and the magnitude of the average state charge 〈q 2 〉 do not change with the RG iterations at all Note that the latter equality is exact (even beyond RG consideration) and equal to the unity 〈q 2 〉 = 1 due to the completeness of the eigenbasis at each ith RG step and for every single realization. Due to Eq. (13) the value of 〈 〉 is completely determined by the zeroth-step cutoff Hamiltonian H 0 or, in other words, by the heads of the wavefunctions.

Equation for effective charges
In order to determine the effective hopping, one can write the equation for χ (R) = 〈|q (R)| 2 〉ν R ( ) = dqP(q, ; R)q 2 6 which is an energy-dependent second q-moment of P(q, ; R), straightforwardly following from Eq. (12) Clearly, the last two delta-functions give after integration −χ ( is the Gammafunction, so we concentrate on the contributions J ± from the first two ones corresponding to δ( − ± ). After changing of integration variables from n and m to w = ( n + m )/2 in both integrals and t ± = ±q n q m f (r) cot(θ /2) in J ± , respectively, one can integrate out the remaining delta-functions and simplify integrands to the identical expressions for J + and J − , The fact that the integration excludes small values of t allows us to simplify the exact relation in the approximation of the small charges. Indeed, assuming an existence of R-dependent cutoff Q(R) for q such that the distribution function P(q, ; R) is exponentially small for q > Q(R) and Q 2 (R) f (R) is small compared to a width of DOS, ν R ( ), for R > R 0 , one can neglect small terms both in the argument of the second P(q, ; R) and in the first brackets getting Here denotes the principle value integration from −∞ to ∞. The function Q 2 (R) in the approximation formulation can be replaced by 〈|q (R)| 2 〉 to obtain a sufficient condition for the validity of the approximation. Thus, the sufficient condition for the approximation to be valid is to have effective hopping t eff (R) = 〈|q (R)| 2 〉 f (R) decaying to zero at infinite distances. As we show below, t eff (R) always behave so in the range of validity of RG scheme, i.e. provided the approximations of effective charges and single resonances are valid. Assuming r.h.s. of Eq. (15) to be sufficiently small even for significantly different R 1 and R 2 , one can replace the finite-difference equation for χ by the following differential one as by taking formally the limit Here, ξ is a natural renormalization scale variable. Solving this equation, one obtains where χ 0 ( ) = χ (ξ = 0) = χ (R = R 0 ), and k 0 ( ) is determined by the expression The condition for replacing the finite-difference equation (15) by the differential one (18) can be rewritten as This is the only validity criterion which explicitly differentiate states by their energy 7 .The condition will later lead to the mobility edge estimation. As seen from (19), the renormalization of χ (ξ) and, thus, of the hopping t i eff , Eq. (8), is non-trivial only if ξ(R) goes to infinity as R → ∞. Otherwise, the effective hopping is proportional to the original one, t eff (R) ∝ f (R). Instead, in the case of non-trivial renormalization χ (ξ) ∼ χ 0 ( )/k 0 ( ) 2 ξ 2 for k 0 ξ 1, and which shows that our original model (2) is dual to the other model, withf (R) ∝ t eff (R) instead of f (R) and localized eigenstates. This result can be applied to any smooth function f (r) and, thus, claims the localization of the bulk spectral states for all long-range ERMs with smooth potential. For example, in a particular case of f (r) = r −a , our result explains the duality of the wave function decay with respect to the critical point a = d observed for d = 1 in [32]. The latter model will be discussed in details in Sec. 3. Note that from Eq. (19) one can easily see that the latter approximation might break down at certain small positive k 0 ( ) due to the presence of a pole which is incompatible with the requirement (21). Thus, at each cutoff R there is a certain energy * determined by the vicinity of the pole k 0 ( * ) = 1/ξ(R) at which χ (ξ) can take unbounded values. It signals that the approximation of single resonances with a given R 0 breaks down for these energies and the states may become delocalized. On the other hand, we know that, due to the general relation 〈q 2 〉 = d χ (R) = 1 working at any R including R 0 the function χ 0 ( ) should be integrable to unity. As a result, χ 0 ( ) should have a sharp peak at the energies * (R 0 ) in the following interval * Indeed, since χ 0 ( ) = ν 0 ( )〈|q (R 0 )| 2 〉, its maximum lies between the absolute maxima of the density of states and the squared effective charge function. The lower bound * min (R 0 ) is of the order of the mean energy which doesn't scale with R 0 , while the upper bound * max (R 0 ) can be estimated as the energy of the trial state with R 0 components all equal to R −d/2 0 with the same sign (zero-momentum plane wave), eventually leading to (24c). This state is chosen as an estimate because it gives the maximal possible value of 〈|q(R 0 )| 2 〉 for the normalized state with R 0 non-zero components Moreover, it is natural to assume that such a state is close to the eigenstate of the model for spatially homogeneous distribution of sites as the fluctuations site positions are averaged out on this zero-momentum plane wave. It is important to note that the presence in such a system of the states with the large effective charge (maybe not only the above mentioned one) which leads to their delocalization causes the localization of the bulk spectral states (similar to the matrix-inversion trick [36]) as the main spectral weight of χ (R) is absorbed by these high-energy states. As we will show in the last section in some models (like in PLE) these delocalized states may be located not at the very spectral edge and this severely questions an alternative cooperative shielding explanation present in the literature [35,52]. Now we are in the position when we are ready to check our approximations, find the range of validity of our method and estimate the size of the wavefunction head R 0 .

Effective-charge approximation
We start with the conditions on the hopping function f (r) required for the effective charge approximation (5). The initial hopping term (4) between states | 0 n 〉 and | 0 m 〉 reads as To go from this form to the approximate one with effective charges we have to assume that f (r i j ) differs only slightly from f (r nm ) for any such i and j that ψ 0 n (i) and ψ 0 m ( j) have significantly non-zero values, i.e. for r ni , r m j < R 0 , Fig. 1. Mathematically, for an isotropic model it gives the following condition: Since the only hopping terms which matter for the RG approach are the ones from the resonant blocks, the distance r nm in the condition should be of the order of a typical distance between counted single resonances. As shown in the next subsection, the next cutoff R i+1 should be chosen to be much smaller than the average distance between single resonances in the full (not truncated) Hamiltonian and, hence, 8 the validity of the effective charge approximation is governed by the condition 9 8 The requirement for the function f (r) to be smooth in points corresponding to the typical distances between resonances, Eq. (28), limits it to have all its sufficiently strong singularities to be located at small distances determined by the zero-cutoff radius, r < R 0 . By the term 'sufficiently strong' we mean such singularities that alter the typical distance between resonances moving it from the value R i+1 towards the vicinity of the pole. Indeed, as soon as the typical resonance is caused by the singular hopping values rather than by R i+1 the corresponding hopping terms cannot be approximated by the effective charges and the whole RG approach fails. 9 The presented condition is sufficient, but far from necessary. An actual necessary and sufficient condition has to deal with meaning of relative fluctuations of f (r) in the R i −vicinity of the typical distance between resonances on the ith step. This fact actually allows the same RG treatment not only for the ERM models with smooth deterministic f (r), but also for the models with hopping terms of the form t i j = (1 + h i j ) f (r i j ) with h i j being a random variable with zero mean and relatively narrow distribution function, f 2 (r i j )〈h 2 i j 〉 r −2d (similar to [37]).
From this relation it is clear that it puts the restrictions not only on the function f (r) but also determines how R i and R i+1 have to be related for a given f (r). For example, in the case of PLE, f (r) ∼ r −a , Eq. (28) gives R i+1 R i . If this restriction contradicts to any other one, then the whole RG approach fails and it may bring the bulk eigenstates of the system to the delocalization.

Single-resonance approximation
Next, we consider the range of validity of the single-resonance approximation for the effective Hamiltonian (5). To justify the approximation, we count the resonances on the ith RG step and estimate the probability for the multiple resonances to occur. For ith RG step, a number of states N i (r) separated by a distance r, R i < r < R, from a certain state with energy and resonant to it can be written as via the probability to have a single resonance at distance r and via the density of states ν i ( ) with energies i n at a given ith RG step. Here ρ(r ) is the average density of sites r in the spatial region of integration. For simplicity we consider the models with uniform spatial density. Thus, we rescale it to unity, ρ(r ) ≡ 1, and omit in further expressions.
Probabilities (30) help us to define the typical probability 10 to have no resonances in the layer R i < r < R, and the typical probability to have exactly one resonance in that layer For the single-resonance approximation to be valid, R i+1 has to be chosen in such a way that the probability to have more than one resonance in the layer, R i < r < R i+1 , is small compared to the probability to have exactly one resonance, i.e.
Assuming the probability p i (r) to be small compared with unity in all points of the layer one can approximately write which finally gives us the smallness, N i (R i+1 ) 1, of the number of resonant states in the layer, Eq. (29), as a requirement. The latter requirement can be written solely via renormalization scales ξ(R i ) and ξ(R i+1 ) in the case of the non-singular R-independent DOS (see Appendix A confirming this assumption) and the decreasing function f (r) Here we used the definition χ (R) = 〈|q (R)| 2 〉ν R ( ), the asymptotic expression for the probability of single resonances p i (r) 2χ (R i ) f (r) for R i , R i+1 1, and the expression (19). Equation (35) together with the condition (21) and Eq. (28) form a complete set of requirements for the RG to be applicable.
The breakdown of the single-resonance approximation occurs when the typical spectrum width is comparable with the effective hopping strength. Indeed, in that case, due to the normalization of the density of states, p i (R i ) ∼ 1 (see Eq. (30)) and, consequently, there is no such R i+1 > R i to satisfy the condition N i (R i+1 ) 1. This result can be intuitively understood as follows: to have a negligible probability of multi-resonance collision we should have low probability of the two-resonance collision as well. So, from this point of view, the first step of our renormalization procedure and an introduction of the zero-cutoff radius were made to remove the singularity of ν R ( ) and made it shallow enough on the scale of the effective hopping.
In the next section we test RG scheme and the above approximations together with wavefunction localization properties for the particular case of the PLE model.

Power-law Euclidean model
To show the validity of our approximations we apply the approach developed in the previous section to the one-dimensional power-law Euclidean (PLE) model and compare our analytical results with numerics. The model is defined by a Hamiltonian (1) with i = 0 As it follows from the previous numerical studies of this model [32], its wave functions show a striking duality, a → 2d − a, of the bulk wave function decay rate. These eigenstates are polynomially localized ψ n,t y p (r m ) 2 ∼ |r n − r m | −2µ with the exponent µ = max{a, 2d − a} for all positive values of the parameter a. Although this model is very similar to the one of Burin and Maksimov (BM) [53] considered also in [36,52,55], the above wavefunction duality in PLE model has not yet been explained theoretically as the matrix inversion trick invented by one of us in [36] breaks down in the absence of diagonal disorder. The RG approach developed in the present paper provides the desired explanation. Indeed, for f (r) ∼ r −a , the non-trivial renormalization occurs for a < d, Eq. (18), giving t eff (r) ∝ r a−2d according to (8), while for a > d renormalization scale ξ converges with R and the standard perturbative approach works giving the polynomial decay of the form r −a . The results for the spatial decay of typical, Eq. (10), and mean, Eq. (11), wave function tails for several cutoff values R corresponding to RG procedure are shown in Fig. 4.
Consider this model in more details. First of all, Eq. (35) in case of PLE model give R i+1 R 2 i , which is compatible with the approximation of effective charges provided R i 1 as R 2 i R i should be valid. It means that the zero-cutoff radius R 0 is finite and large compared to the unity. Fig. 4 provides the following estimate of the cutoff radius R 0 ∼ 30, which is in full agreement with the above consideration. However, our RG approach is actually applicable only for a > 0: for negative values of a the approximation of single resonances breaks down since the original hopping f (r) is no longer a decreasing function. Nevertheless, according  10) and (11) (written in panels as equations). The right column shows that the validity of RG scheme for typical wavefunction decay can be extended also to some spatially increasing (though unphysical) hopping.
to numerical results, this breakdown of the RG approach doesn't lead to the delocalization or even to the aforementioned duality breaking for typical wavefunction spatial decay. This fact may be caused by the destructive interference of the resonances or, in other words, by the higher-order corrections to the perturbation series. After determining the validity range of RG, we check numerically the fact about the density of states stated in the Appendix A. According to the RG approach, the function ν R ( ) barely depends on the cutoff radius R, a < d: (37) and it converges to a non-singular function. Both these statements are clearly seen from Fig. 5 supporting the analysis done in the Appendix A.
Next, although the estimates and numerical results presented above in this section justify the RG approach for PLE model in general, they do not provide any information about the approximations we made going from Eq. (15) to Eq. (18), i.e., the small-charge approximation and the approximation of the finite-difference equation by the differential one. To check that the effective charges are indeed behave according to Eq. (19), we calculate it explicitly for a set of different cutoff values, see Fig. 6. The lower row of panels shows the function 〈|q (R)| 2 〉ξ 2 (R) which does not depend on the cutoff value and collapses to a universal curve with good accuracy. Moreover, the insets show that the above collapse works relatively well until the maximum of 〈|q (R)| 2 〉, but not only in the bulk of the spectrum. Finally, our theory predicts that, for a < d, 〈|q (R)| 2 〉 as a function of must have a sharp maximum at * ma x ∼ R d−a 0 with the magnitude of the order of R d 0 , see (24) and the corresponding discussion. By combining these two estimates we get the one describing the energy dependence of the maximal magnitude with increasing cutoff R As shown in the insets to Fig. 7, this is exactly the case, at least for d = 1.  25) and (38). The insets show the same maximal points of 〈|q (R)| 2 〉 in linear scale with the power-law black dashed fitting curves coinciding with the ones in the main panels. The data is averaged over 10 3 realizations.
One may notice that the maximum of the effective charge 〈|q | 2 〉 in Fig. 7 occurs at the edge of the spectrum for a < d/(d +1) = 1/2, while for a > 1/2 the maximal energy corresponds to the constant R-independent value of 〈|q | 2 〉 = O(1) 11 .The explanation of this is based on the fact that on top of the trial delocalized state there are rare states localized at few adjacent sites. For R d 0 sites the minimal distance between such sites, typical for each disorder realization, is given by r min = R −d 0 12 and thus the energy of the state localized on this pair of states scales as ∼ f (r min ) ∼ R d a 0 and at a > d/(d + 1) these states will form the edge of the spectrum. The corresponding effective charges for these states are given by the expression lim →∞ 〈|q | 2 〉 = 2.
As mentioned in the first section, the presence of the delocalized states with large energies scaling with the system size (or cutoff value) causes the localization in long-range Euclidean matrices in the similar way as in the models with diagonal disorder and translation-invariant hopping terms [36] due to the leakage of most of the charge spectral weight to large energies (and measure zero of states). The lesson which one should take from this is the following: it is not the ground (or anti-ground) state with the energy diverging with the system size which matters for the localization of the bulk spectrum (like in the cooperative shielding approach [35,52]), but the presence of high-energy delocalized states (with high effective charges) do this job. The latter high energy states do not need to be at the very spectral edge, see the right panel of Fig. 7. 11 Note that the case a < d/(d + 1) = 1/2 characterized by the delocalized eigenstates at the very spectral edge is an artefact of finite statistics in our numerical simulations. Indeed, in the renormalization group written for the infinite system with a certain cutoff, see, e.g., the bottom left of Fig. 3, has to be determined by the infinite number of states localized at few very close sites, which, in turn, form the very spectral edge. In numerics instead of the cutoff R i of the infinite matrix we diagonalize full R i × R i matrices removing many two-site localized states. Another effect of such numerics is that it forces the localization radii to be not larger than R i at each ith RG step and reduces the corresponding finite-size effects for the wavefunction tails (shown in Fig. 1). 12 This estimate is given by solution of the equation R d 0 P(r min ) = 1, with the distribution of distances between adjacent sites, homogeneously distributed in d-dimensional space with unit density, given by Poisson formula, P(r) ∼ r e −r .

Conclusion and discussions
To sum up, in this work we develop a generic renormalization group (RG) approach applicable to a wide range of Euclidean random matrix models, which shows localization of the bulk mid-spectrum states and provides the wave function decay for these states in an explicit form, Eq. (8). The range of validity of the above statements is governed by three conditions: the applicability of effective charge approximation (28) restricting the hopping potential f (r) to be smooth, the single-resonance approximation (35) which is satisfied, e.g., for the bounded monotonically decaying function f (r), and the slow probability density evolution condition (21) which is deeply interconnected with the single resonances approximation.
The above mentioned requirements allow us to get rid of any small parameter, which is crucial for the standard RG approach [33,34,53,54], and show the renormalization to the localization for all spectral bulk eigenstates. This localization is solely caused by the drastic spectral flow of the renormalized effective hopping to high-energy delocalized states (forming measure zero of all states in low dimensions d ≤ 2). The developed RG has many similarities to the so-called matrix inversion trick developed by one of us in [36] and complements and extends it to the case of Euclidean matrices with off-diagonal disorder (with or even without diagonal disordered part).
Moreover, the developed approach shows the equivalence between Euclidean models and translation invariant models with diagonal disorder and smooth hopping potential f (r) not only in the localization properties, but also in the spatial decay of bulk mid-spectrum eigenstates. We believe that this equivalence can be generalized to non-smooth and even anisotropic hopping which is recently under the spotlight [52,55], but this is the topic of further investigation. structure of the distribution function P(q, ; R) and expand it up to at least the first order of the Taylor series assuming q 2 n q 2 m f 2 (r)/t to be small. As a result, we get the differential equation After substituting here the approximate expression for χ from Eq. (19), we find that, in case of large ξ(R) for large R, the derivative ∂ R ν R ( ) ∼ R d−1 f 2 (R)/ξ 3 (R), i.e., As long as f (R) goes to zero with increasing R, the r.h.s of the latter expression does the same. This concludes that, for R R 0 , the function ν R ( ) is saturated and only slightly differs from its limiting value ν( ) = ν ∞ ( ).