Localization and fractality in disordered Russian Doll model

Motivated by the interplay of Bethe-Ansatz integrability and localization in the Richardson model of superconductivity, we consider a time-reversal symmetry breaking deformation of this model, known as the Russian Doll Model (RDM), and implement diagonal on-site disorder. The localization and ergodicity-breaking properties of the single-particle spectrum are analyzed using a large-energy renormalization group (RG) over the momentum-space spectrum. Based on the above RG, we derive an effective Hamiltonian of the model, discover a fractal phase of non-ergodic delocalized states -- with the fractal dimension different from the paradigmatic Rosenzweig-Porter model -- and explain it in terms of the developed RG equations and the matrix-inversion trick.

In this paper, motivated by the interplay of the BA integrability, localization, and level repulsion in the Richardson model, we consider the Russian Doll model bringing TRS breaking into the game, along with diagonal disorder. As in previous localization studies of the Richardson model, we focus on the single-particle sector of the RDM, which has much in common with the many-body physics including the tower of high-energy ground state solutions. We also consider the generalization of RDM in terms of scaling of the coupling constant. In the original RDM, the coupling scales as N −1 , while we consider more general scaling N −γ/2 in analogy with the Rosenzweig-Porter model [12]. The latter model also consists of all-to-all hopping terms, but the couplings are given by i.i.d. Gaussian random variables with the standard deviation N −γ/2 .
The Rosenzweig-Porter model is known to host an entire phase of non-ergodic (socalled fractal) eigenstates in the range 1 < γ < 2, squeezed between the ergodic (γ < 1) and Anderson localized (γ > 2) phases [13]. The non-ergodic phase is characterized by the only energy scale Γ, large compared to the level spacing δ ∼ 1/[N ρ(E)] and small compared to the bandwidth ∼ 1/ρ(E) of the spectrum, where ρ(E) is the density of states. This energy scale is given by the standard Fermi's Golden rule and it determines the fractal dimension 0 < D = 2 − γ < 1 of the wave-function support set. Later, several other models with similar fractal [14][15][16][17][18] and multifractal [19][20][21][22][23][24] phases have been suggested in the literature. In all these models, it has been shown that the wavefunction structure is determined mostly by the diagonal elements, while the hopping terms provide a certain Breit-Wigner level broadening Γ. 1 . For the Richardson model, the standard Fermi's Golden rule result fails to describe the localization properties correctly due to the presence of strong correlations between the coupling of different sites. In the case of localization -which survives for any coupling strength even beyond the convergence of the locator expansion (like in the Richardson model [4,5] and other long-range fully correlated models [26,27]) -one can use a so-called matrix-inversion trick [15] or develop a strong-disorder spatial RG [28,29].
For RDM, in this work we show that increasing the coupling does lead to the delocalization of most of the eigenstates, and therefore, both the above methods that work only in the localized phase, are not applicable. At the same time, the standard Fermi's Golden rule approximation (1) fails due to the strongly correlated coupling terms. Therefore, our goal here is to develop another analytical method to describe localization and ergodicitybreaking properties of RDM. We base our approach on the RG flow, similar in spirit to the one used for disorder-free RDM for γ = 2 [7], but generalize it to the momentum space. In order to double-check the RG approximations, we also generalize the above-mentioned matrix-inversion trick to the case of any unbound spectrum of the disorder-free coupling term j mn , and show that the effective Hamiltonian obtained by this method is statistically equivalent to the one calculated from the RG.
By going back to the coordinate basis, we derive the effective Hamiltonian with significantly reduced correlations, which is readily tractable with the Fermi's Golden rule approximation (1). The effective Hamiltonian makes it possible to elaborate the localization properties of the single-particle states. Using a combination of the effective Hamiltonian and Fermi's Golden rule, we find that single-particle eigenstates in the disordered RDM demonstrate fractal properties, emerging at the same Anderson localization point γ = 2 as the Rosenzweig-Porter model. However, the non-ergodic phase is prolonged to smaller values of γ, i.e. γ = 0, and the corresponding fractal dimension D, which we determine exactly analytically, also deviates from the one in the Rosenzweig-Porter model and equals to D = 1 − γ/2.
The remainder of the paper is organized as follows. In Sec. 2 we explicitly describe the disordered Russian Doll model. Next, in Sec. 3 we calculate the spectrum of the disorder-free RDM, describe it in terms of energy stratification [17], and calculate the localization properties of the energy-stratified states in the momentum basis. In Sec. 4 we derive an effective Hamiltonian representation for RDM using a high-energy RG in the momentum space. Section 5 represents the generalization of the matrix-inversion trick introduced in [15] in order to make it applicable to the description of the delocalized states and confirm its equivalence to the above RG by comparing the results for the effective Hamiltonian. In Sec. 6 we provide the analytical results leading from the structure of the effective Hamiltonian, supported with numerical simulations. The conclusion and outlook are given in Sec. 7.

Model
In this work, we focus on the single-particle sector of the Russian Doll model with on-site disorder ε n and generalized coupling amplitude j mn ∼ N −γ/2 . The single-particle Hamiltonian in the coordinate basis is the N × N random matrix where 1 ≤ m, n ≤ N . Here, the generalized coupling j mn scales as power −γ/2 of the system size N , as opposed to N −1 , and the on-site potentials ε n are given by Gaussian i.i.d. variables ⟨ε n ⟩ = 0, The above-mentioned symmetric coupling g and the TRS breaking parameter h are parameterized by the angular variable θ as follows For simplicity, we consider periodic boundary conditions and define the distance d(m, n) between sites m and n with a sign: if the shortest route from m to n is clockwise (counterclockwise), the distance is positive (negative), see Fig. 1, This allows us to determine an effective magnetic flux θ, threading the loop m − n − m and equal for each link between any pair of sites m and n.
Note that according to Anderson resonance counting [13,15,30,31], a general principle of Anderson localization in long-range models, a measure one subset of the states in this model are localized for γ > 2, irrespective of any correlations or TRS breaking. Consequently, all the properties present in the Richardson model or Rosenzweig-Porter model at γ > 2 -such as the Lorentzian power-law profile of the eigenstates versus ε i (sometimes called frozen multifractality) and the power-law Chalker scaling of the wavefunction overlap, absent in the short-range Anderson models -are also present in the disordered Russian Doll model. Therefore, unless mentioned otherwise, we focus on the range 0 < γ < 2 in the further sections. (2)-(4). Different colors of vertices stand for the disorder potential ε n , while the coloring of the edges from the topmost vertex demonstrate different phases of hopping terms with the same amplitude: red color stands for e iθ , blue -for e −iθ , and black dashed line corresponds to the real hopping 1.

Energy stratification of the spectrum and ergodicity-breaking in momentum space
In this section, we focus on the spectrum of the disorder-free RDM. Using the property of its stratification, we analyze the localization and the ergodicity-breaking properties of large-energy eigenstates of the corresponding disordered RDM. Indeed, the coupling/hopping matrix j is translation invariant j mn = j m−n and can be thus diagonalized in the basis of plane waves with the spectrum indexed by an integer |p| ≤ N/2 From this spectrum one can immediately see that • For the Richardson model (θ = 0), the spectrum is (N −1)-fold degenerate, E p̸ =0 = 0, with the only non-zero energy level E 0 ∼ N 1−γ/2 . It is this level which is responsible for the localization of the other N − 1 eigenstates orthogonal to it in the disordered Richardson model for γ < 2 [4,5].
• Even in the general case of θ ̸ = 0, the levels with non-zero and even p = 2k remain small, |E 2k | < N −γ/2 , for odd N , and are zero for even N . Later, we focus on the case with even N to neglect the small amplitude of these levels. However, the levels with odd p = 2k + 1, for any finite θ, are as significant as E 0 . For |k| ≪ N one can write E 2k+1 ∼ sin θ 2N 1−γ/2 π(2k + 1) . (8) Note that the transition between Richardson model and RDM (see the above cases) occurs for the TRS breaking parameter θ c ∼ W N −(1−γ/2) . For θ < θ c , the largest of the energy levels E 2k+1 , namely E 1 , becomes smaller than the diagonal disorder amplitude W and thus becomes hybridized with the rest of the zero modes by the disorder. As N → ∞, this transition occurs at θ c → 0. Therefore, the Richardson model is an exceptional point, owing to the discontinuity in the behaviour of RDM as θ → 0 and the Richardson model at θ ≡ 0 in the thermodynamic limit N → ∞ 2 . The transition at another special point θ = π/2 is continuous as the only level E 0 = 0 goes to zero at that value.
In the many-body sector of the Richardson model and RDM, there is one BCS-like ground state, or a whole hierarchy of such states, and the gap(s) from these states to the rest becomes extensive at γ < 2. The single-particle sector of these models demonstrates the same structure of gapped or energy stratified levels, and moreover, the number of such levels also scales similarly with the system size N . In the Richardson model (and other long-range fully correlated models [26][27][28]) as well as in RDM, the energy stratified levels are special: they form a measure zero subset of all the spectral states, but give the main contribution to the hopping term j mn . The most high-energetic of these states are barely affected by the disorder term, and thus, stay non-ergodic in the momentum basis due to their extensive diagonal energy. This leads to both, the ergodicity of these states in the real space, and the fact that they give the main contribution to the hopping term.
Indeed, the disorder term ε n (3) in the momentum-space basis (6) plays the role of scattering between plane waves (or hopping) with translation-invariant Gaussian i.i.d.
amplitudes J p−q = 1 N n e 2πin(p−q) N ε n , with zero mean and the variance scaled down with the system size Thus, the corresponding representation of RDM in the momentum-space is a translationinvariant realization of the Rosenzweig-Porter ensemble with a special choice E p of the diagonal disorder. For this model introduced in [15] it is known that the Fermi's Golden rule is applicable and gives the following broadening of the Breit-Wigner approximation for the eigenstate (see, e.g., [19,32,33]) Here, C is an unimportant normalization constant and we labelled the high-energy eigenstates with disorder-free energy E p assuming smallness of the broadening Γ p with respect to it. One should note that, unlike the Rosenzweig-Porter model, the RDM in the momentum space has a highly inhomogeneous density of states (DOS) ρ(E p ), p = 2k + 1, given by where we have taken into account that the disorder ε n ∼ W hybridizes the levels as soon as the disorder-free version of DOS |dp/dE p | goes above its bare disorder counterpart |dn/dε n | ∼ 1/W . The support set ∆p occupied by the eigenstate (11) in the momentum space can be found using the condition which implies non-ergodic behavior as soon as ∆p ∝ N D(p) scales as a fractional D(p) < 1 power of N . As soon as |E p+∆p − E p | ≃ |E p | (or ∆p ≃ p), the condition (13), using (8), (10), and (12), leads to the number p * of energy-stratified states which are non-ergodic in the momentum basis: The energies of these states are barely affected by disorder as E p are extensive, E p ≫ W for |p| ≪ p * . Thus, our criterion is consistent with the so-called Mott's principle (see, e.g., [15]), which claims that as soon as the bare diagonal energy E p of a state is large compared to the spectral width W of the hopping term J p−q , this state is non-ergodic in the corresponding (momentum) basis. Note that for θ < θ c ∼ N −(1−γ/2) only the state p = 0 is non-ergodic (localized) in the momentum space. Note also that the support set ∆p is limited from above by p < p * ∼ N 1−γ/2 . Thus, from the condition D(p) < 1, we see that Eq. (14) is valid until p * ≪ N , i.e. for γ > 0. Henceforth, we will mostly focus our considerations to this parameter interval, 0 < γ < 2.

Large energy RG in the momentum space & Effective Hamiltonian
In the paper [7], the authors consider renormalization over the matrix size N in the disorder-free RDM with linearly increasing diagonal terms ε n ∼ n. Each RG step involves the removal of one row and one column corresponding to the largest diagonal element ε N . In particular, the RG considered in [7] can be described as follows: 1. Start with the matrix of size N 0 and reduce its size by one at each step.
2. For this, at each step, take the largest absolute diagonal energy (ε N or ε 1 ), and assuming it to be large with respect to the rest of the levels and the hopping terms, resolve the eigenproblem with respect to the site i = N corresponding to the level ε N : where m ̸ = N and j mn (1) is calculated by one RG step with j mn (0) = j mn .
3. Next, assume ε N + j N N − E ∼ W and using the ratio W/δ = N one obtains cyclic RG equations.
In the disordered RDM, the procedure discussed above fails as the maximal diagonal matrix-element does not correspond to the maximal (or minimal) index, which breaks down the self-similar structure of the matrix at further RG steps (see, e.g., [34]). However, one can consider a renormalization group analogous to Eq. (17) by taking into account the large diagonal terms in the momentum basis, where the diagonal energies E p (7) and the hopping terms J p−q (9) satisfy the inequality (15), since The corresponding equation for the hopping term at rth step of the renormalization, for removal of the level with momentum p r , is given by while the diagonal terms stay the same Of course, the described renormalization works only when Eq. (18) is valid, i.e. at least for γ < 3, which is trivially satisfied in our interval of the interest, 0 < γ < 2. Further, the removal of the level with the largest |E p | proceeds in the following order, for s up to s = r ≤ N/4 Here we should warn the reader that both θ = 0 and θ = π/2 have been considered slightly differently. In the vicinity of the Richardson model, θ ≲ N −(1−γ/2) , the only level satisfying Eq. (18) is E 0 , thus we can consider only one renormalization step r = 1. On the contrary, in the vicinity of θ = π/2, the level E 0 invalidates (18), thus we should start with s = 1. However, as we will see in Eqs. (27) and (31), the latter choice does not change the results. Further, we plan to find the optimal number of RG steps needed for writing the effective Hamiltonian for the bulk spectral states, E ∼ O(1), with suppressed correlations. The effects of the energy-stratified states will be taken into account by the RG. In order to find the effective Hamiltonian, in subsection 4.1, we first simplify the RG flow (19) focusing on the leading contributions by order of magnitude. Next, in subsection 4.2, we rewrite the Hamiltonian in the coordinate basis in order to find the optimal number r of RG steps needed to minimize the broadening Γ, found by Fermi's Golden rule from the effective Hamiltonian.

Simplification of RG (19)
In order to simplify the RG equation (19), here we show that its main contribution is given byJ where Indeed, this takes into account the renormalization (19) itself, but neglects the renormalization of the hopping terms J p,q (r) in the sum S p,q (r). As we show in Appendix A, for r smaller than the above approximation works well, leading to |S p,q (2r)| ≪ |J p−q |, and the difference between J p,q andJ p,q is at most of the order |S p,q |. Note that the value r * * , corresponding to momentum p * * = 2r * * − 1 ≫ p * according to Eq. (21), is large compared to the number p * of energy-stratified levels (14) for γ > 0. Thus, one can take into account all the high-energy states within the above RG flow.

Effective model in the coordinate basis.
Now we are in a position to rewrite the effective renormalized model (20) and (22) in the coordinate basis in order to estimate the fractal dimension of the eigenstates. For this purpose, we separate our renormalized Hamiltonian into four terms where a p,q,l = J p−2l+1 J 2l−1−q − J p+2l−1 J −2l+1−q and p, q ̸ = p s , with 0 ≤ s ≤ r, and p s are from Eq. (21). The discrete Fourier transform of the above terms takes the form where we have replaced the summation over p, q ̸ = {p s } by complementary sums over the whole interval and over p s in either or both variables. The first summation is given simply by the initial (not truncated) Fourier transform. After some straightforward algebra and neglecting subleading corrections (both given in Appendix B), one obtains the following renormalized Hamiltonian in the coordinate basis at 2rth step of the RG flow, with 1 ≤ r ≤ N/4 and an unimportant constant c 5 Generalization of the matrix-inversion trick for Russian Doll model Here we present an alternative way to derive the effective Hamiltonian of RDM in the momentum space, analogous to Eq. (25), which is free from the approximations of the above RG. For this purpose, we generalize the matrix-inversion trick developed in [15]. The main idea of the matrix-inversion trick is as follows: given a hopping matrix with large eigenvalues E p , one adds to it, the identity matrix multiplied by a certain constant E 0 , and inverts this matrix as follows In this way, the large energies E p of the hopping matrix, that provide the dominant contribution to the hopping, can be sent to the denominator without changing the basis.
Thus, the effective model can be treated with perturbation theory as soon as the parameter E 0 is chosen so as to avoid any resonances In models with one-sided unbounded growth of the spectrum E p (like the Burin-Maksimov model [15] where E p<p * ≫ 1 are large and positive), one can avoid having singularities in the denominator, E p +E 0 , by choosing E 0 < − min p E p ∼ O(1), and demonstrate wave-function localization. However, in the RDM, the spectrum is unbounded on both sides (see Eq. (8) for positive and negative k ≪ N ) and does not have finite gaps at finite energies in the thermodynamic limit. Thus, one cannot find a suitable E 0 ∼ O(1) to avoid the divergence arising from the inversion of E p + E 0 terms.
In order to obtain convergent terms while applying the matrix-inversion trick, one has to invert only a part of the spectrum given by high-energy levels E p<pr where, for simplicity, we choose E 0 = −E and use |E |p|<pr | ≫ E.
Observe that Eq. (29) corresponds term-by-term with Eq. (25). Indeed, • the first term is equivalent to H (1) after neglecting subleading terms like i mn in (65) of Appendix B, • the part of the second term with p = p 0 = 0 corresponds to EH (2) mn /ε m after neglecting g m,0 r/N subleading terms in (25), • the rest part of the second term with p = p s ̸ = 0 corresponds to EH (3) mn /ε m after neglecting g m,0 r/N subleading terms in (25), • while the last term is just equal to H (4) .
To sum up, this result shows that all the approximations performed in the previous section in order to derive the effective Hamiltonian, Eq. (27), either lead to subleading corrections or to prefactors E/ε m ∼ O(1) of order unity. As the matrix-inversion trick provides merely a different representation of the exact eigenproblem without any approximations, the equivalence between Eqs. (25) and (29) confirms the applicability of the effective Hamiltonian (27) in the whole range of parameters of interest, 0 < γ < 2.

Results
Now we are ready to calculate the non-ergodic properties of eigenstates based on effective Hamiltonian (27) and Fermi's Golden rule approximation (1). Similar to the matrix-inversion trick [15], the Fermi's Golden rule for each effective Hamiltonian with a certain r gives an upper bound for the fractal dimension via the broadening Γ n (2r) (see definition 30). Therfore, in order to find the true fractal dimension, one should make the upper bound strict by finding the optimal r = r opt that minimizes Γ n (2r). Below, we implement the optimization procedure analytically, and later verify our result for the fractal dimension using numerical calculations of eigenstate statistics.

Analytical results -Optimization of fractal dimension
Similar to the Rosenzweig-Porter model, we expect the model given by the effective Hamiltonian (27) to exhibit only fractal (defined via level broadening in (33)) and not multifractal states (where multifractal dimensions D q are parameterized by the order q of the wave-function moment, see Sec. 6.2). Therefore, in order to calculate the fractal dimension D of eigenstates in RDM, we use Fermi's Golden rule analogous to (1), but with the effective hopping term from the renormalized Hamiltonian (27) Γ n (2r) = 2π Taking the energies and DOS in the bulk of the spectrum, E n ∼ ε n ∼ W and ρ(E n ) ∼ 1/W , one obtains, for each term of the effective Hamiltonian, the following expression 3 The first term (corresponding to H (3) with p = p 0 ) is subleading for all |θ| ≫ N −(1−γ/2) and r ≫ 1. Formally, in the vicinity of θ = π/2 this term diverges, but as we discussed earlier, the inequality (18) should be satisfied in order to write this term. For similar reasons, in the vicinity of θ = 0, the divergence of the third term ∼ 1/ sin θ can be ignored.
After neglecting the first term for finite θ, the rest of the three terms give the optimal value of r corresponding to the minimal level broadening with certain constants c and c ′ of order one. Note that this optimal value r opt satisfies the condition (24), r ≪ r * ∼ N 1−γ/3 , for all γ > 0. The validity condition |J p,q (r) −J p,q (r)| ≪ |S p,q (r)| for the above used RG, leading to r ≪ N (3−γ)/4 from Eq. (56) in Appendix A, is satisfied for γ > 1. However, even for r ≳ N (3−γ)/4 , corresponding to 0 < γ < 1, J p,q (r) −J p,q (r) gives at most the same-order contribution as S p,q (r) and affects only the numerical prefactors c and c ′ .
In the bulk of the spectrum, as the mean level spacing is δ = 1/(ρ(E)N ) ∼ W/N , the fractal dimension for the typical wave function in this case should be given by the ratio which is different from the one in the Rosenzweig-Porter model [13]. This result is also confirmed by numerical calculations below, where we define the fractal dimension via the inverse participation ratio (but not as the number of sites where the wave function has significantly non-zero values within the Breit-Wigner approximation). Note that in cases where the Fermi's Golden rule does not work for the initial problem, the number p * of the energy-stratified levels, Eq. (14), determines the fractal dimension via the expression p * ∼ N D (see [17] for more details). The effective renormalized Hamiltonian (27) in this case is given by where W/Γ ≃ N/r opt ∼ (W/ sin θ)N γ/2 and we neglect the subleading terms for simplicity. Note that the Hamiltonian (34) is equivalent to the one of power-law random banded matrices [35] with a bandwidth b = W/Γ and diagonal disorder W rescaled as ∼ N γ/2 .
• if none of b and W scales with N , the system hosts genuinely multifractal states, with D q ∼ b at b ≪ 1 and 1 − D q ∼ 1/b at b ≫ 1 [35], • the scaling b ∼ N γ/2 , W = O(1) sends the system into the ergodic phase, with D = 1, • the scaling W ∼ N γ/2 , b = O(1) leads to the localization, D = 0, • whereas the case of RDM, corresponding to the simultaneous scaling of both parameters W ∼ b ∼ N γ/2 , gives fractal eigenstates described by the Fermi's Golden rule (30) 4 .
The bulk eigenstates of such an effective Hamiltonian should be given by two contributions: • First, due to the presence of Rosenzweig-Porter-like long-range hopping terms at |m − n| ≪ W/Γ, the wave-function should have a contribution of a Lorentzian profile versus ε n with the width |E m − ε n | ∼ Γ [19,32,33] |ψ • Second, similar to the power-law banded random matrices |j mn | 2 ∼ |m − n| −2a at the critical power a = 1, there should be the multifractal proliferation of the wavefunction maxima given by the resonances. However, unlike the power-law banded random matrix case, in the RDM, the N -scaling of the cutoff W/Γ ∼ N γ/2 , at which this power-law comes into play, significantly reduces the number of resonances: and does not affect the fractal dimension of the system, given by (see e.g., Eqs. (30)(31) in [36])

Numerical results -Spectrum of fractal dimensions, level statistics, and wave-function decay
In numerics, we proceed free of approximations and consider exact diagonalization of the initial model (2), calculating eigenstates ψ En (m) in the coordinate basis and eigenvalues E n . We analyze the data using probes of multifractality while focusing on the mid-spectrum states. For this purpose, we consider two relevant measures of eigenfunction statistics based on the distribution of amplitudes P (|ψ E (n)| 2 ).
First, we look at the spectrum of fractal dimensions, defined as the power f (α) of the scaling of the distribution, The spectrum f (α) is a kind of large deviation function showing the tails of the distribution of ln |ψ E (n)| 2 far from its typical (most probable) value ln |ψ E (n)| 2 = −α 0 ln N .
It has a bunch of properties (normalization condition of the probability distribution f (α) ≤ 1, f (α 0 ) = 1, or the wave-function normalization f (α) ≤ α, f (α 1 ) = α 1 ), among which we specifically mention the symmetry [37], originally discovered in [38], relating peaks (small α < 1) and tails (large α > 1) of the wave-function. This symmetry is known to work for non-ergodic extended states, fractal or multifractal. In particular, for the Rosenzweig-Porter model, as shown in [13], f (α) takes a simple linear form for γ ≥ 1, with an additional point f (0) = 0 for γ > 2: The above spectrum satisfies the symmetry (40) for γ < 2. The second widely used probe of multifractality that we consider is the inverse participation ratio I q defined via the moments of eigenstates, as follows, where the q-dependent exponents D q are called fractal dimensions. In the ergodic phase, D q = 1 for all q, whereas in the case of localization, D q = 0 for q > 0. The intermediate values, 0 < D q < 1, correspond to non-ergodic extended states, which can be either multifractal, where D q is represented by a strictly decaying function of q, or fractal, where D q = D does not depend on q, at least for q > 1/2. The relation between D q and f (α) is given by the saddle-point approximation for disorder-averaged moments, and the Legendre transform: where the last equality is the definition of α q . From the above definition, one can immediately see that α 1 , determining the wavefunction normalization condition f (α 1 ) = α 1 , gives the limiting value D q→1 where the latter equality is given by the symmetry (40). Here α 0 determines the most typical wave-function amplitudes (39) and consequently, gives the maxima, f (α 0 ) = 1.
Using the above relation (45) and our prediction for value of the fractal dimension, Eq. (33), one can find the spectrum of fractal dimension in the RDM, similar to (41), Alternatively, one can derive the above expression from the Breit-Wigner approximation (11) using the broadening (33), see e.g. [33]. Numerically, as one cannot achieve infinite system sizes, both the spectrum of fractal dimensions Eq. (38) and the fractal dimensions themselves Eq. (42), are calculated for different finite system sizes and extrapolated to the thermodynamic limit N → ∞ (see Appendix C). The spectrum of fractal dimension at finite system sizes can be extracted directly from the histogram over α (see, e.g., [13,15,39,40]), while the inverse participation ratio is given simply by Eq. (42).   Figure 2 shows the spectrum of fractal dimensions f (α) extracted from the numerical simulations. One can see that the extrapolation procedure gives the correct normalization value f (α 0 ) = 1, and good agreement with the analytical prediction (46).
One can use either the inverse participation ratio (42), or its relation to f (α) (45), in order to extract the fractal dimension, see Fig. 3. Agreement amongst these two measures and with the analytical prediction (33) confirms our analytical derivation and the symmetry (40) of f (α). Slight deviations from the prediction in the localized phase are caused by finite-size effects, which are enhanced close to the Anderson transition.
Note that, according to the analytical prediction, the fractal dimension is not homogeneous across the spectrum and this is also confirmed numerically, see Fig. 4: The fractal dimension in the spectral bulk (more than 90 % of states at the considered system sizes) shows the above fractal behavior for 0 < γ < 2, while at the edges, the high-energy states become ergodic with D 2 → 1.
In addition to the previous two measures, we also consider the eigenvalue statistics using a so-called adjacent level gap ratio (defined in [41,42]) and the wave-function spatial decay (first used in [27]). The ratio statistics considers Localization corresponds to Poisson level statistics, characterized by the absence of level repulsion, and leads to mean r = ⟨r n ⟩ = r P = 2 ln 2 − 1 ≃ 0.3863. The ergodic randommatrix prediction corresponds to the Wigner surmise [6] and is given by r GOE ≃ 0.5307, for orthogonal symmetry and r GU E ≃ 0.5996 for the unitary one. For the fractal phase of the Rosenzweig-Porter kind, the gap ratio still shows the Wigner-Dyson value in the entire delocalized phase, γ < 2, while non-ergodicity is visible only at higher order gap ratios (see, Figure 4: The fractal dimension D 2 (E n ) versus energy index n in RDM model for γ = 0.75, 1, 1.25, θ = 0.25, and N = 2 10 , 2 12 , 2 14 . (insets) the same data for D 2 (E n ), shown vs energy E n and zoomed close to the right spectral edge. The data both for D 2 (E n ) and E n are averaged over 1000 realizations of disorder for each eigenvalue index separately.
e.g., [29,43,44]) or in other measures like the spectral form factor, level compressibility [13] or the power spectrum [45][46][47]. Figure 5 shows the conventional ratio statistics (47) across the spectrum in disordered RDM for γ = 0.75, 1, 1.25. One can see that, similar to the Rosenzweig-Porter model, the disorder-averaged ratios in the bulk of the spectrum follow the Wigner-Dyson unitary value r GU E , while the spectral edges states correspond to special r values, first going up to the equidistant spectral average r = 1, and then to the nearly degenerate levels r ≃ 0.
Our last measure, the wave-function spatial decay [15,27,29,48], uncovers the structure of the eigenstates predicted by the effective Hamiltonian (34) and Eq. (35). In order to observe the wave-function spatial decay, we plot it versus the diagonal energy differences |ε m − ε n |, provided these energies are ascendingly ordered ε n < ε n+1 , see Fig. 6.
Indeed, according to (35), for sites n with the diagonal energy ε n close to the eigenvalue E m , |ε n − E m | ≲ Γ, the wave-function has a fractal structure |ψ Em (n)| 2 ∼ N −D like in the Rosenzweig-Porter model. This is confirmed by the Lorentzian form of the wave-function decay in Fig. 6, similar to the Rosenzweig-Porter results [48].

Conclusion and Outlook
In summary, we consider the localization and ergodicity-breaking properties of eigenstates in the disordered Russian Doll model, with generalized amplitude of the coupling strength ∼ N −γ/2 . In this regard, we develop RG flow based on the renormalization of high-energy states in the momentum basis, and provide an effective Hamiltonian, valid for any number of renormalized high-energy states within the considered parameter range.
In addition, we validate our result for the effective Hamiltonian and confirm that the approximations we used are valid to leading order in the parameter r/r * * , i.e. the number of the RG steps normalized by the sub-extensive upper bound r * * ∼ N 1−γ/3 , and all subsequent corrections are subleading. We do so by generalizing the matrix-inversion trick to cases where the hopping-spectrum is dense at all finite energies and grows without bound on both sides of the spectrum. Unlike the localized case, described by the matrix-inversion, the above spectral structure corresponds to the phases of delocalized eigenstates that can be understood using the above generalization of the matrix-inversion trick (cf. [17]).
Based on effective Hamiltonian, we find the fractal dimension of eigenstates in the bulk of the spectrum. We see that RDM has a delocalised non-ergodic phase for an extended range of parameter values (0 < γ < 2) compared to the Rosenzweig-Porter model (1 < γ < 2), and with a different fractal dimension D = 1 − γ/2.
Note that, as the fractal properties in RDM barely depend on the time-reversal symmetry breaking parameter θ ̸ = 0, the Richardson model provides an example of an exceptional point, since the limiting behavior θ → 0 of the Russian Doll model does not correspond to that of the Richardson model, θ = 0. Unlike the Burin-Maksimov model [15,26,27] -where the symmetry-breaking parameter destroys the BA integrability, and thus, one expects to see discontinuous behaviour -this work demonstrates that broken time-reversal symmetry in RDM, which preserves the BA integrability for any finite N , may lead to a similar discontinuity.
It would be interesting to look for non-ergodicity in the ensemble of twisted XXX models with random inhomogeneities, using their relation to the RDM model. Another interesting direction for further study concerns the derivation of the generalized RDM, that is, Richardson model with broken time-reversal symmetry and the generalized hopping term scaling, and the investigation of its fractal properties. It is an important direction of research as the generalized Richardson model describes the superconducting phase of the mixture of Sachdev-Ye-Kitaev and Fermi-Hubbard models [49]. Moreover, we show that the Russian Doll model provides an example of a Bethe-Ansatz integrable model with delocalized non-ergodic eigenstates already in the single-particle sector. This raises concerns about the relation between Bethe-Ansatz integrability and localization and opens a new avenue in this research direction. A Estimation of smallness of a parameter S p,q (2r) and J p,q (r)− J p,q (r) Within the condition (24), r ≪ r * * = N 1−γ/3 , and the sum S p,q (r), Eq. (23), is small compared to J p−q where a p,q,l = J p−2l+1 J 2l−1−q −J p+2l−1 J −2l+1−q due to (9) has zero mean and the following variance Using this, the above sum of random variables can be approximated via its variance (as it have the zero mean) Finally, this gives the following estimate for S p,q (2r) at r ≫ 1 due to γ < 3 and (24). In a similar way one can estimate the difference l p,q (r) = J p,q (r) −J p,q (r) , l p,q (1) = 0 .
Indeed, using Eqs. (19), (22), and (48) one immediately obtains where we neglected the quadratic term (S p,p k (r) + l p,p k (r)) (S p k ,q (r) + l p k ,q (r)) due to its smallness. Further we estimate by the order of magnitude the parameter l by rewriting the above equation for l in the continuous form and neglecting the difference between variables with different indices dl(r) dr ≃ J E pr (S(r) + l(r)) .
Solving this ordinary differential equation in the variable x one obtains For x ≪ 1 one can immediately see that |l(r)| ∼ xS(r) ≪ S(r).
In the opposite limit of x ≳ 1 one can only bound |l(r)| ≤ S(r) using the condition for y ′ = (x − y) ≥ 0 in the integrand In this case one cannot neglect l(r) with respect to S(r), but can absorb it to S(r) if the numerical prefactors are not important. Thus, in the main text we consider both above cases of x ≪ 1 and x ≳ 1.
B Derivation of the effective Hamiltonian (25) in the coordinate basis (27) In Eq. (25) we separate our renormalized Hamiltonian in four terms where we introduce the notation a p,q,l = J p−2l+1 J 2l−1−q − J p+2l−1 J −2l+1−q and p, q ̸ = p s , with 0 ≤ s ≤ r, and p s are from Eq. (21). The discrete Fourier transform of the above terms takes the form where we replaced the summation over p, q ̸ = {p s } by the complemented sums over the whole interval and over p s in either or both variables. The first summation is given just by the initial (not truncated) Fourier transform.
Here we will calculate all these terms one-by-one. The first term written in the above four sums takes the form where the first term corresponds to the diagonal disorder, the second one is given by while the last two terms are symmetric with respect to each other by the Hermitian conjugation with p ′ = q and q ′ = p and J −p = J * p . Let's calculate, first, I mn by shifting the summation over p to k = p − q The second sum i mn in Eq. (61) can be found after the same shift (65) and estimating the following sum with the random phase approximation and the central limit theorem for r ≫ 1 leading to random variable g m of the order of one. After neglecting of the small terms g m r/N with respect to ε m we obtain for the first term In the last equality we approximate the sine factors of the last term by linear functions when their arguments are small compared to one. The second term H (2) m,n splits in the product of two equivalent sums, giving in the same approximation as for i mn The third term H m,n within the same approximation reads as In the last case we assumed that sine of the large argument is more or less equivalent to (−1) l . In both latter derived equations we again used (66) and neglected these terms with respect to ε m . The last term is given by the truncated initial hopping term which we do not split into the above four sums (60): Again using the same approximation for sine of the large argument as (−1) k we will obtain The first bracket in the first case corresponds to the small argument of the sine, k ≤ N/[4π(m − n)], while the rest terms correspond to the large sine argument.
To sum up, in the coordinate basis at 2rth step we have the following estimate for the renormalized Hamiltonian given in Eq. (27) H m,n (2r) ∼ ε m δ mn + ε m ε n N 2−γ/2 cos θ + +   C Extrapolation of the multifractality measures to the thermodynamic limit N → ∞ Here we remind the standard extrapolation procedure for the spectrum of fractal dimensions (see, e.g., [13,15,16,27,39]) and for the fractal dimensions D q [37].
For the first one we express the multifractal spectrum f (α, N ) at finite system size N with certain constants A (k) α depending on α. The latter expression can be derived using the definition Eq. (38) and extracted directly from the histogram over α [13,15,39,40].
Here and further we stick to the simplest linear in 1/ ln N behavior, which is typical for the models with fractal eigenstates [13,15,17]. The corresponding finite-size f (α, N ) and extrapolated f (α) curves are given in Fig. 7 for 50 % of mid-spectrum states. As an additional marker of the extrapolation quality we check that the normalization condition, max α f (α) = f (α 0 ) = 1, of the probability distribution P(α) is satisfied.
The finite-size fractal dimension is defined by the formula (42) D q (N ) = ln I q /(1 − q) ln N , with the generalized inverse participation ratio (IPR), In order to avoid the parasitic contributions from measure zero of special eigenstates we focus on the typical averaging of the IPR both over the disorder and the eigenstates I q,typ = e ⟨ln Iq⟩ ∼ N −(q−1)Dq,typ and omit the subscript "typ" for brevity.  Fig. 7. Different symbols in the extrapolation of α q correspond to different percentage P of the deviation from the maximum of the function f (α q ) − qα q used for the extrapolation.
As the main contribution to IPR is given by the scaling exponent D q and the prefactor c q similarly to (73) one obtains The extrapolation of D q (N ) vs 1/ ln N extracted from I 2 and from α 0 and α 1 is shown in Fig. 8. Here in order to diminish finite α-bin size for extracting α q we fit f (α) − qα close to its maximum with a parabolic fit and associate α q with the maximal position of this fit. The fitting interval, α − ≤ α ≤ α + , is determined by the deviation from the maximal value f (α q ) − qα q to f (α ± ) − qα ± = (f (α q ) − qα q ) P , with the percentage P = 0.9, 0.93, 0.95, 0.97 shown in the legends of Fig. 8. One can straightforwardly see that 7 % change in P affects the extrapolation by at most the same amount. The same procedure is done for α 0 (N ) and α 1 (N ) mentioned in Eqs. (39) and (44).