# The Wishart–Rosenzweig–Porter random matrix ensemble

Victor Delapalme, Leticia F Cugliandolo, Grégory Schehr, Marco Tarzia, Davide Venturelli

► To cite this version:

Victor Delapalme, Leticia F Cugliandolo, Grégory Schehr, Marco Tarzia, Davide Venturelli. The Wishart–Rosenzweig–Porter random matrix ensemble. Journal of Statistical Mechanics: Theory and Experiment, 2026, 2026 (1), pp.013301. 10.1088/1742-5468/ae29f6. hal-05462596

HAL Id: hal-05462596

https://cnrs.hal.science/hal-05462596v1

Submitted on 16 Jan 2026

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
PAPER • OPEN ACCESS

# The Wishart–Rosenzweig–Porter random matrix ensemble

To cite this article: Victor Delapalme et al J. Stat. Mech. (2026) 013301

View the [article online](#) for updates and enhancements.

You may also like

- Overlap gap and computational thresholds in the square wave perceptron
  Marco Benedetti, Andrej Bogdanov, Enrico M Malatesta et al.

- Deterministic simplicial complexes
  S N Dorogovtsev and P L Krapivsky

- The correlation time of a super-Poisson renewal process characterized by waiting time distributions with the locally thickened tail
  Aleksey V Kargovsky, Olga A Chichigina and Alexander A Dubkov
PAPER: Disordered systems, classical and quantum

# The Wishart–Rosenzweig–Porter
random matrix ensemble

Victor Delapalme<fnref n="1" />,<fnref n="2" /> ID, Leticia F Cugliandolo<fnref n="3" />,<fnref n="4" /> ID,
Grégory Schehr<fnref n="3" /> ID, Marco Tarzia<fnref n="1" />,<fnref n="4" /> ID
and Davide Venturelli<fnref n="1" />,* ID

E-mail: davide.venturelli@sorbonne-universite.fr

Received 23 October 2025

Accepted for publication 21 November 2025

Published 16 January 2026

<figure><img src="image_2.png" /></figure>

CrossMark

Online at stacks.iop.org/JSTAT/2026/013301

https://doi.org/10.1088/1742-5468/ae29f6

**Abstract.** In recent years, the Rosenzweig–Porter (RP) ensemble, obtained by adding a diagonal matrix with independent and identically distributed elements to a Gaussian random matrix, has been widely used as a minimal model for the emergence of fractal eigenstates in complex many-body systems. A key open question concerns the robustness of its phase diagram when the assumption of independent and uncorrelated entries is relaxed—an assumption that simplifies its analysis but is generally violated in realistic quantum systems. In this work, we take a first step in this direction by considering a deformed Wishart (rather than Gaussian) random matrix, which we dub the ‘Wishart–RP’ ensemble. Using perturbation theory, as well as the cavity and replica methods and the Dyson Brownian motion approach, we characterize its phase diagram and localization properties. Remarkably, we show that the level compressibility, which quantifies

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹ Laboratoire de Physique Théorique de la Matière Condensée, CNRS-UMR 7600, Sorbonne Université, 4 Place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote n="2">² Department of Mathematics, King's College London, Strand, London WC2R 2LS, United Kingdom</footnote>

<footnote n="3">³ Laboratoire de Physique Théorique et Hautes Energies, CNRS-UMR 7589, Sorbonne Université, 4 Place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote n="4">⁴ Institut Universitaire de France, 1 rue Descartes, 75005 Paris, France</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
spectral correlations in the fractal phase, coincides with that of the Gaussian RP model, thereby extending the universality conjectured in Venturelli *et al* (2023 *SciPost Phys.* **14** 110) beyond the fully uncorrelated setting. We confirm our results with numerical tests.

**Keywords:** cavity and replica method, random matrix theory and extensions, Anderson transition and Anderson localization, ergodicity breaking

Contents

<table><tr><td>1. Introduction</td><td>3</td></tr><tr><td>2. The WRP ensemble</td><td>5</td></tr><tr><td>2.1. Summary of the main results</td><td>6</td></tr><tr><td>3. Perturbation theory: Mott's criteria for localization and ergodicity</td><td>10</td></tr><tr><td>3.1. Spectrum of fractal dimensions</td><td>12</td></tr><tr><td>4. The cavity method</td><td>15</td></tr><tr><td>4.1. Density of states, local density of states and phase diagram</td><td>16</td></tr><tr><td>4.2. Density–density correlations and level compressibility in the intermediate regime</td><td>19</td></tr><tr><td>5. The replica method</td><td>22</td></tr><tr><td>5.1. Density of states</td><td>22</td></tr><tr><td>5.2. Full counting statistics</td><td>23</td></tr><tr><td>5.3. Level compressibility</td><td>27</td></tr><tr><td>6. Dyson Brownian motion</td><td>28</td></tr><tr><td>7. Conclusions</td><td>31</td></tr><tr><td>Acknowledgments</td><td>32</td></tr><tr><td>Appendix A. Scaling function of the level compressibility</td><td>32</td></tr><tr><td>Appendix B. Details of the replica calculation</td><td>33</td></tr><tr><td>B.1. Density of states</td><td>33</td></tr><tr><td>B.2. Full counting statistics and level compressibility</td><td>37</td></tr><tr><td>Appendix C. Details of the Dyson Brownian motion calculation</td><td>39</td></tr><tr><td>C.1. First-order perturbation theory for the eigenvalues</td><td>39</td></tr><tr><td>C.2. Second-order perturbation theory for the eigenvalues</td><td>40</td></tr><tr><td>C.3. First-order perturbation theory for the eigenvectors</td><td>41</td></tr><tr><td>C.4. Second-order perturbation theory for the eigenvectors</td><td>42</td></tr><tr><td>References</td><td>43</td></tr></table>
# 1. Introduction

Random matrix theory (RMT) [1–3] has long been an invaluable framework for describing and understanding complex physical systems. Its power lies in universality: many RMT results are largely independent of the specific distribution of matrix entries, making them applicable across a broad variety of physical contexts. One prominent example is quantum chaos and its breakdown. In particular, RMT underpins our current understanding of quantum ergodicity, formalized through the eigenstate thermalization hypothesis (ETH) [4, 5].

In recent years, a large body of work has suggested that disordered interacting quantum systems may violate ETH under certain conditions. A paradigmatic example is many-body localization ([6, 7]; see [8–13] for recent reviews), which occurs at strong disorder. More generally, the emergence of multifractal eigenstates—i.e. states that do not uniformly explore the accessible Hilbert space, thereby violating ETH and quantum ergodicity, hence often called non-ergodic<fnref n="5" />—has been identified as a robust feature of the phase diagram of such systems [14–30], at the origin of several of their unconventional properties. To capture these phenomena, simple random-matrix models have been proposed and extensively studied [31–69]. The philosophy behind these approaches is that RMT offers analytically tractable models capable of explaining general and universal features of systems that otherwise resist analytic treatment.

The most prominent such model is the (generalized) Gaussian Rosenzweig–Porter (GRP) ensemble [31, 70]. It is defined as the sum of two independent $N \times N$ matrices, namely a diagonal random matrix $A$ with independent and identically distributed (i.i.d.) entries, and a matrix $B$ from the Gaussian orthogonal ensemble (GOE) with random elements with zero mean and variance of $O(1)$,

$$ \mathbf{H}_{\text{GRP}} = \mathbf{A} + \nu N^{-\gamma/2} \mathbf{B}. \qquad (1) $$

In the prefactor $\nu N^{-\gamma/2}$, the parameter $\nu$ is of $O(1)$. The physical interpretation is intuitive: each site of the reference space (matrix index) corresponds to a configuration of the system with a random on-site energy drawn from $A$, while transitions between configurations are mediated by Gaussian-distributed amplitudes from $B$. The spectral properties of the model are controlled by the parameter $\gamma$. This model provides a prototypical example of a system that exhibits an intermediate non-ergodic extended phase ($1 < \gamma < 2$), characterized by fractal eigenstates and unconventional spectral properties, lying between a fully delocalized phase ($\gamma < 1$) and a fully Anderson-localized phase ($\gamma > 2$) [31].

For this reason, the RP model and its generalizations have recently received renewed attention as a playground to explore the nature and properties of non-ergodic extended states [31–67, 71–73]. A key open question is the extent to which the spectral properties of the RP ensemble are robust under modifications of the distribution of the Hamiltonian

<footnote n="5">⁵ Similarly, fully-delocalized eigenstates that satisfy ETH are often called ergodic [4, 5]. Strictly speaking, however, in the context of single-particle non-interacting problems such as the one considered here, the concept of ergodicity is not sharply defined (especially because the eigenvalues lack the extensive scaling typical of interacting systems). Nevertheless, following the common usage in the literature, throughout this paper we will use the term non-ergodic as a synonym of (multi)fractal.</footnote>
matrix elements. One of the main motivations of this work is to address precisely this issue. In particular, we will focus on the spectral compressibility, which is closely related to the two-point spectral correlation function. This quantity, the definition of which is recalled in section 2.1, displays distinct behaviors in the delocalized regime (with level repulsion described by RMT) and in the localized regime (with uncorrelated Poisson statistics).

In a recent work [61], some of us derived the exact scaling function describing the crossover between these two regimes in the intermediate phase of the GRP model. There, we showed that this scaling function is universal with respect to the distribution of the entries of **A**. More recently, a generalization of the RP model, where the matrix **B** is drawn from the Lévy ensemble [42, 43, 68] (with i.i.d. entries and power-law tails), was also shown<fnref n="6" /> to yield the same universal scaling function [69].

In this paper, we introduce another generalization of the RP model: the Wishart–RP (WRP) ensemble, in which **B** is drawn from the Wishart ensemble [74], a fundamental class of random matrices [1, 2]. An ensuing difference with respect to the standard GRP case is that in the Wishart case the matrix **B** is positive definite, with all positive eigenvalues, thereby producing an asymmetric rightward shift of the energy levels of **A**. More importantly, the entries of **B** are statistically dependent, featuring non-zero higher-order correlations, even though the pairwise correlations vanish (see below). Using a combination of complementary approaches, including perturbation theory, the cavity method, the replica method and Dyson Brownian motion, we obtain the full phase diagram of the WRP ensemble and study its spectral properties in its different phases. We then compute the spectral compressibility and the two-point correlation function in the intermediate regime, showing that the same universal crossover scaling function as in the GRP model emerges, although the elements of **B** are not independent. These predictions are fully supported by exact diagonalization numerics.

Our results support the idea that, at least for models where the intermediate phase is characterized by a simple fractal (rather than multifractal) spectrum with compact mini-bands (see section 2.1 for a precise definition), the crossover from RMT universality to Poisson statistics is genuinely universal. This calls for a systematic numerical analysis of the crossover function in more realistic many-body systems.

Finally, we stress that the implications of our results extend beyond toy models of ergodicity breaking in disordered quantum systems. Indeed, variants of the WRP model are relevant in other scientific contexts as well. For example, in denoising problems, whose aim is to recover a signal hidden in a noisy covariance matrix of many correlated time series, this ensemble provides a mathematically controlled framework to characterize the statistical structure of noise at different scales, and to separate it from the meaningful signal [75–79]. Other examples are generative machine learning models: from the perspective of generative diffusion models, which gradually transform an initial random state through a learned reverse diffusion process, the WRP model naturally emerges in the late stages of the backward diffusion process (see e.g. [80]).

<footnote n="6">⁶ Strictly speaking, the authors of [69] actually computed the two-point density–density correlation function $\langle\rho(\omega_1)\rho(\omega_2)\rangle_c$, which is related to the level compressibility via the relation in equation (58) below.</footnote>
The rest of the presentation is organized as follows. In section 2 we define the WRP ensemble, and in section 2.1 we present a summary of our results. Next, we employ multiple methods, namely the perturbation theory in section 3, the cavity method in section 4, the replica method in section 5, and the Dyson Brownian motion in section 6, to derive the phase diagram, the average density of eigenvalues and the spectral compressibility. We support our results with exact diagonalization data, which we discuss throughout the paper. Finally, in section 7, we conclude and note some directions for future research.

## 2. The WRP ensemble

The WRP ensemble is defined as the sum of two $N \times N$ random matrices $\mathbf{A}$ and $\mathbf{B}$:

$$ \mathbf{H} = \mathbf{A} + \nu \mathbf{B}, \qquad (2) $$

where $\nu$ is a constant of $\mathcal{O}(1)$. Here, $\mathbf{A}$ is a random diagonal matrix $A_{ij} = a_i \delta_{ij}$, where the $a_i$'s are $N$ independent and identically distributed random variables with probability density $p_a$ (we assume $\langle a \rangle = 0$, and $p_a(0) > 0$). In equation (2), $\mathbf{B}$ is a Wishart matrix:

$$ \mathbf{B} = M^{-\gamma} \mathbf{W} \mathbf{W}^T, \qquad (3) $$

where $\mathbf{W}$ is an $N \times M$ rectangular matrix with i.i.d. Gaussian distributed entries, i.e.

$$ \rho_W[\mathbf{W}] = \frac{1}{(2\pi)^{NM/2}} \prod_{i=1}^{N} \prod_{\ell=1}^{M} e^{-W_{i\ell}^2/2}, \qquad (4) $$

such that $\langle W_{i\ell} \rangle = 0$ and $\langle W_{i\ell}^2 \rangle = 1$. We define $c = N/M \le 1$, and we restrict ourselves to the case in which $N$ and $M$ are of the same order, i.e. $c$ is of $\mathcal{O}(1)$. The elements of the matrix $\mathbf{B}$ are given by

$$ B_{ij} = M^{-\gamma} \sum_{\ell=1}^{M} W_{i\ell} W_{j\ell}; \qquad (5) $$

hence, one has

$$ \langle B_{ii} \rangle = M^{1-\gamma}, \quad \langle B_{ij} \rangle = 0, \quad \langle B_{ij}^2 \rangle = M^{1-2\gamma} \quad \text{for } i \neq j. \qquad (6) $$

Note that the matrix $\mathbf{B}$ in equation (5) is positive definite, i.e. with all positive eigenvalues. In the large-$N$ limit, the eigenvalues of $\mathbf{B}$ are distributed according to the celebrated Marčenko–Pastur law [1, 2]. Since $\langle \text{Tr} \mathbf{B} \rangle \propto NM^{1-\gamma} \propto N^{2-\gamma}$, it follows that the eigenvalues of $\mathbf{B}$ (and, consequently, the support of the Marčenko–Pastur distribution) scale as $N^{1-\gamma}$. In particular, they grow with $N$ for $\gamma < 1$, decrease with $N$ for $\gamma > 1$, and remain of $\mathcal{O}(1)$ for $\gamma = 1$.

At first sight, the ensemble (2) appears similar to the model introduced in [50], where the matrix $\mathbf{B}$ is likewise constructed as a sum of independent projectors. However, there
<figure><img src="image_4.png" /><figcaption>Figure 1. Phase diagram of the WRP model (see section 2.1 for its description).</figcaption></figure>

is an important distinction between the two cases: in [50], the matrix **B** is defined so that its energy levels are independent random variables, and therefore follow Poisson statistics. In contrast, in our setting, **B** has a rotationally invariant statistical weight, and its energy levels exhibit the universal correlations of RMT at all scales. Although the phase diagrams of the two models share qualitative similarities (see section 3.1 below), this distinction results in fundamentally different spectral correlations. In this sense, the WRP model studied here is more closely related to the standard GRP ensemble than to the model of [50].

Nonetheless, the WRP ensemble also differs from the standard GRP ensemble in a key aspect: the entries of **B** are statistically dependent. In particular, from equation (5) it is straightforward to show that, albeit the connected two-point correlations vanish, the correlations between *n*-tuples of matrix elements with repeated (but different) indices are given by

$$ \langle B_{i_1 i_2} B_{i_2 i_3} \cdots B_{i_n i_1} \rangle = M^{1-n\gamma}. \qquad (7) $$

Note that, since $M = N/c$, the parameter *c* has the effect of tuning these correlations (at fixed *N*), the maximally correlated case corresponding to the $c \to 0$ limit.

## 2.1. Summary of the main results

Here, we present a brief summary of the main findings of our work. The first concerns the phase diagram of the model, shown in figure 1, which features three distinct phases separated by two transition points, as detailed below:

* $\gamma < 1$: In this regime and for large *N*, the matrix **A** is subleading with respect to **B**, and thus the latter completely dominates the spectral properties of **H**. The average density of states (DoS) is given by the Marčenko–Pastur distribution [81], with eigenvalues and support growing as $N^{1-\gamma}$, up to subleading finite-*N* corrections, which depend on **A**. This corresponds to a fully delocalized phase with Wigner–Dyson statistics.

* $\gamma > 3/2$: Here, the matrix **A** dominates over **B** when *N* is large, and fully controls the spectral properties of **H**. The eigenvalues of **H** are close to the diagonal entries of **A**, up to small perturbative corrections induced by **B**, which are not strong enough to hybridize more than $O(1)$ energy levels. This corresponds to a fully localized phase, in which the eigenstates are localized around those of **A**, and the energy levels obey Poisson statistics.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Illustration of the different scales in the spectrum in the intermediate phase of the WRP model ($1 < \gamma < 3/2$). The blue window illustrates the concept of mini-bands in the spectrum, extending over the scale of the Thouless energy $E_T$.</figcaption></figure>

* $1 < \gamma < 3/2$: This intermediate regime is the most interesting one. The matrix **B** can still be treated as a perturbation, but in this case it hybridizes energy levels on a scale that is parametrically much larger than the spectral gap (i.e. the average level spacing), yet still much smaller than the total bandwidth. As a result, the average DoS is still given by the distribution of the diagonal entries of **A**, but correlations between energy levels emerge on an energy scale known as the *Thouless energy* $E_T$, involving a number of hybridized states that grow with $N$ as $N^{2-2\gamma}$ (see figure 2). This corresponds to a partially ergodic *fractal* phase, analogous to that of the standard GRP model. The eigenvectors are partially delocalized on a fractal support set that grows with $N$ as $N^D$, with $D = 3-2\gamma$, but is much smaller than $N$.

The subleading corrections to the average DoS (which converges to the Marčenko-Pastur distribution for $\gamma < 1$ and to $p_a$ for $\gamma > 1$) are given by the Zee formula [82, 83] for the average density of eigenvalues of the sum of two random matrices, for which we provide a derivation within the replica approach (see equations (49), (71) and (72) below).

The second and main result of this work concerns the (super)universal behavior of the *level compressibility* (or, equivalently, of the density-density correlation function) in the intermediate phase on the scale of the Thouless energy. The level compressibility $\chi(E)$ is a simple indicator of the degree of level repulsion and is defined as follows [84]. Let

$$ \rho^{(e)}(\lambda) = \frac{1}{N} \sum_{i=1}^{N} \delta(\lambda - \lambda_i) \qquad (8) $$

denote the ‘empirical’ eigenvalue density, and

$$ I_N[\omega_1, \omega_2] = N \int_{\omega_1}^{\omega_2} d\lambda \rho^{(e)}(\lambda) \qquad (9) $$

denote the number of eigenvalues $\lambda_i$ lying in the interval $[\omega_1, \omega_2] \subseteq \mathbb{R}$, which is a random variable. Denoting $E = (\omega_2 - \omega_1)/2$ the width of the energy window, and $E_0 = (\omega_2 + \omega_1)/2$ its middle point, the level compressibility is defined as
The Wishart–Rosenzweig–Porter random matrix ensemble

$$ \chi(E) \equiv \frac{\kappa_2(E)}{\kappa_1(E)} = \frac{\langle I_N^2 \rangle - \langle I_N \rangle^2}{\langle I_N \rangle} = \frac{\langle I_N^2 \rangle_c}{\langle I_N \rangle}, \qquad (10) $$

where $\kappa_1$ and $\kappa_2$ are the first two cumulants of $I_N$. For Poisson statistics$^7$, one has $\kappa_2(E) \simeq \kappa_1(E)$, and then $\chi(E) \simeq 1$. On the other hand, for a rigid spectrum like that of the Wishart matrix $\mathbf{B}$, the mean number of eigenvalues behaves for small $E$ as $\langle I_N \rangle \propto N\rho(E_0)E$, where

$$ \rho(\lambda) \equiv \langle \rho^{(e)}(\lambda) \rangle, \qquad (12) $$

while $\langle I_N^2 \rangle_c \propto \ln(N\rho(E_0)E)$ [84]. Hence, in this case, one finds $\chi(E) \to 0$ for $E \gg [N\rho(E_0)]^{-1}$ (but still much smaller than $E \sim \mathcal{O}(1)$).

As we show below, in the intermediate fractal phase of the WRP ensemble, the level compressibility assumes a scaling form. This function, which we plot in figure 3, describes the crossover from the universal RMT behavior at small energy separations to Poisson statistics at large separations, and it reads

$$ \chi \left( y = \frac{E}{E_T} \right) = \frac{1}{\pi y} \left[ 2y \operatorname{atan}(y) - \ln(1 + y^2) \right]. \qquad (13) $$

The small- and large-$y$ asymptotics of this scaling function are given in equation (64) below. Remarkably, this expression is identical to that of the GRP ensemble, both when $\mathbf{B}$ is real and symmetric (as in the GOE), and when it is Hermitian (as in the Gaussian unitary ensemble (GUE)).

In [61], some of us showed that, on the scale of the Thouless energy, the level compressibility of the GRP model is insensitive to the specific distribution of the i.i.d. diagonal entries of $\mathbf{A}$ (up to a rescaling by $E_T$ and provided $p_a(0) > 0$). More recently, the same crossover function was found in the Lévy–RP ensemble, where $\mathbf{B}$ is a Lévy matrix with entries drawn from a power-law distribution [69] (see, however, footnote 6). Here, we demonstrate (using both the cavity method and the replica approach) that the very same scaling function (13) also governs the crossover in the WRP ensemble, which differs in two key aspects from the Gaussian and Lévy RP ensembles: (i) $\mathbf{B}$ has positive-definite eigenvalues, inducing an asymmetric rightward shift of the spectrum of $\mathbf{A}$, and (ii) its matrix elements are not independent (see equation (7)).

We have tested this theoretical prediction by performing exact numerical diagonalization of large random matrices from the WRP ensemble around the Thouless energy $E_T$, varying the matrix size $N$ while keeping the ratio $c = N/M$ fixed, in the intermediate phase ($\gamma = 1.25$). In figure 3, we compare the numerical results with the analytical scaling prediction (13), finding an excellent agreement. In figure 4(a) we show that the

$$ \chi(E) = 1 - \frac{\langle I_N[\omega_1, \omega_2] \rangle}{N} = 1 - \int_{\omega_1}^{\omega_2} d\lambda \rho(\lambda). \qquad (11) $$

We thus generically expect $\chi(E) \sim 1$ for small $E$, and $\chi(E) \to 0$ for large $E$.

<footnote>$^7$ However, for energy separations of the order of the total bandwidth, $\chi(E)$ actually decreases from 1 to 0 at very large $E$ (see e.g. appendix A of [61]). In fact, one has $\kappa_1 = \langle I_N[\omega_1, \omega_2] \rangle$ and $\kappa_2 = \langle I_N[\omega_1, \omega_2] \rangle (1 - \langle I_N[\omega_1, \omega_2] \rangle / N)$, and thus the level compressibility reads</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Level compressibility of the WRP ensemble. The numerical data obtained via exact diagonalization (green symbols) are compared to analytical predictions in the distinct regimes of the model, as detailed below. First, at small energy scales (comparable to the mean level spacing $\Delta = 1/(N\rho)$), the level compressibility decreases in accordance with the universal behavior $\chi_{GOE}$ of Wigner–Dyson statistics (see equation (91)). In this regime, the spectrum is said to be rigid because the variance of the number of eigenvalues in an interval is small due to level repulsion (and thus, so is the level compressibility). Conversely, for intervals of $O(1)$, the level compressibility matches that of i.i.d. random variables (11), meaning that there are no correlations between energy levels at this scale. Finally, for intervals close to the Thouless energy $E_T$, the numerical curve approaches the (super)universal crossover function $\chi_T$, which connects the two regimes; see equation (13). Note that the replica method provides an analytical prediction (yellow symbols; see section 5) that is valid for all regimes, except the one at $E \lesssim \Delta$, where it breaks down. For this plot, we used $\gamma = 1.25$, $c \approx 0.998$, $\nu = 1$, $N = 2000$, and $p_a$ is the standard Cauchy distribution.</figcaption></figure>

range of validity of (13) increases with the system size when the energy separation is measured in units of the Thouless energy. The crossover function (13) is also found to be independent of the parameter $c$ that controls the correlation strength; see figure 4(b).

These observations are quite remarkable since the same crossover scaling function is found in many distinct random matrix ensembles, for which indeed the level compressibilities do not necessarily coincide for $E \gg E_T$ or $E \ll E_T$. This finding strongly suggests that such crossover functions are in fact fully and genuinely universal, and that such universality originates from the structural properties of the model, rather than from the specific choice of matrices **A** and **B**, at least across all models in which the intermediate phase is fractal rather than multifractal. It also motivates numerical studies of this crossover function in more realistic many-body disordered quantum systems that exhibit (multi)fractal phases, to test whether the same universal behavior emerges in that context as well.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Numerical data for the level compressibility of the WRP model at the Thouless energy scale, obtained by exact diagonalization of large random matrices, for (a) different values of *N* with *c* = 0.98 fixed, and (b) different values of *c* with *N* = 20000 fixed. We chose *γ* = 1.25 and *ν* = 1; moreover, we took *p*<sub>*a*</sub> to be uniform between −1 and 1, which allowed us to average the level compressibility over multiple non-overlapping intervals (which are statistically equivalent since the density of states is constant). The dashed line indicates the universal scaling form χ<sub>T</sub> given in equation (13).</figcaption></figure>

**3. Perturbation theory: Mott's criteria for localization and ergodicity**

The simplest and most intuitive way to analyze the phase diagram of the model is
through first- and second-order perturbation theory for the eigenvalues and the eigen-
vectors of **H** in equation (2), which yield the so-called *Mott's criteria* for localization
and ergodicity [85]. These have been successfully applied to the GRP model and its
generalizations (see e.g. [45] for a detailed explanation).

In the limit where the off-diagonal matrix **B** is absent, all eigenvectors |ψᵢ⟩ are
trivially localized on a single site, i.e. |ψᵢ⟩ = |i⟩ (where |i⟩ is the position basis), with
corresponding eigenvalues λᵢ = aᵢ. The first criterion, known as *Mott's criterion for
localization*, states that Anderson localization around a single matrix index occurs when
the average level spacing Δ = 1/(Nρ) is much larger than the tunneling amplitude
between different indices. The second criterion, known as *Mott's criterion for ergodicity*,
provides a sufficient condition for ergodicity. The idea is to estimate the average escape
rate Γ of a particle localized on a given site using Fermi's golden rule, and to compare it
to the spread of energy levels. When the average spreading width Γ is much larger than
the energy bandwidth, the different indices are fully hybridized: starting from a given
site, the wave packet spreads to any other site with the same energy on a timescale of
O(1).

In equations, Mott’s criterion for localization reads: ⟨|*H*<sub>*ij*</sub>|⟩ ≪ Δ, where Δ = 1/(Nρ) is the average gap, ρ being the average spectral density at the considered value of energy. Since for γ > 1 the eigenvalues of the matrix **B** are of the order N<sup>1-γ</sup> and they all vanish in the large-N limit, the average DoS is asymptotically given by the probability distribution of the entries of the matrix **A**, namely ρ(λ) = p<sub>*a*</sub>(λ). For γ < 1, instead, the entries of **A** are much smaller than the eigenvalues of **B** in the large-N limit, and the average DoS is given by a Marčenko–Pastur distribution with support
growing as $N^{1-\gamma}$. We thus obtain that the average level spacing is

$$ \begin{cases} \Delta = 1/(Np_a(\lambda)) = 1/(Mcp_a(\lambda)), & \text{for } \gamma > 1, \\ \Delta \propto N^{-\gamma}, & \text{for } \gamma < 1. \end{cases} \qquad (14) $$

Since the spectral properties of the model do not depend on the specific energy at which they are probed, in the following, without loss of generality, we focus on the center of the energy band, i.e.

$$ \begin{cases} E_0 = 0 + \mathcal{O}(N^{1-\gamma}), & \text{for } \gamma > 1, \\ E_0 \propto N^{1-\gamma}, & \text{for } \gamma < 1. \end{cases} \qquad (15) $$

From the definition of the model, one immediately obtains

$$ \langle |H_{ij}| \rangle = \nu M^{-\gamma} \left\langle \left| \sum_{\ell=1}^{M} W_{i\ell} W_{j\ell} \right| \right\rangle = \nu \sqrt{\frac{2}{\pi}} M^{1/2-\gamma}. \qquad (16) $$

Mott's criterion then states that Anderson localization close to the eigenvectors of A occurs provided that

$$ \langle |H_{ij}| \rangle = \nu \sqrt{\frac{2}{\pi}} M^{1/2-\gamma} \ll \frac{1}{cp_a(0)M} \quad \Rightarrow \quad \gamma > 3/2. \qquad (17) $$

Applying second-order perturbation theory to the eigenvalues and the Fermi golden rule, we now compute the average bandwidth $\Gamma$ that corresponds to the energy window within which hybridization occurs:

$$ \Gamma = 2\pi p_a(0) \left\langle \sum_j H_{ij}^2 \right\rangle = 2\pi p_a(0) N \nu^2 M^{1-2\gamma} = 2\pi p_a(0) \nu^2 c M^{2-2\gamma}. \qquad (18) $$

The quantity $\hbar\Gamma$ can be interpreted as the bandwidth that can be reached in a time of $\mathcal{O}(1)$ from a given site $i$, and is often called the Thouless energy $E_T$. This implies that the eigenvectors within this energy window are hybridized by the Wishart perturbation.

For $1 < \gamma < 3/2$, such energy band decreases with the system size as $E_T \propto N^{2-2\gamma}$, but is still much larger than the mean level spacing $\Delta \propto N^{-1}$. This entails that the system is not Anderson localized; nevertheless, $E_T$ remains much smaller than the total bandwidth, which is of $\mathcal{O}(1)$. This signifies that the particle can only explore a subextensive portion of the total Hilbert space.

The Anderson localization transition occurs when $E_T$ becomes smaller than the mean level spacing, i.e. for $\gamma > 3/2$. This implies that the average escape time from site $i$, defined as $\Delta t \equiv \hbar/E_T$, grows at least linearly with $N$, and thus the eigenfunctions remain localized on $\mathcal{O}(1)$ sites. In contrast, the transition to the fully delocalized phase takes place when $E_T$ becomes of the order of the total bandwidth, i.e. for $\gamma \le 1$. For $\gamma < 1$, one has that $E_T$ is much larger than the total bandwidth $N^{1-\gamma}$. Hence, starting at site $i$, a wave packet can reach any other site in a time of $\mathcal{O}(1)$, corresponding to full delocalization.
In the intermediate phase $1 < \gamma < 3/2$, the support set of the eigenvectors (i.e. the number of sites hybridized by the perturbation) is given by the ratio between the width of the hybridized energy window and the average gap between adjacent energy levels. It therefore scales as

$$ \frac{E_T}{\Delta} \sim N^D, \quad D = 3 - 2\gamma. \qquad (19) $$

The partially extended but fractal eigenstates are thus linear combinations of $N^D$ localized states associated with nearby energy levels. The wave function of one of these eigenstates can thus be represented as [35, 42, 61]

$$ |\psi\rangle \approx \sum_{i:|a_i|<E_T} c_i N^{-D/2} |i\rangle, \qquad (20) $$

with $c_i$ being a Gaussian random variable with zero mean and variance 1. These eigenstates give rise to the so-called *mini-bands* in the local spectrum (see figure 2). The Thouless energy thus corresponds to the energy window within which Wigner–Dyson-like spectral correlations (and in particular level repulsion) are established.

All the moments $I_q$ of the wave-function coefficients (the so-called generalized inverse participation ratios) behave as

$$ I_q = \sum_i |\langle i|\psi\rangle|^2 q \propto N^{D_q(1-q)}, \qquad (21) $$

which defines the fractal dimensions $D_q$. The schematic picture presented in equation (20) corresponds to a situation in which all fractal dimensions $D_q$ are degenerate and equal to $D$. In section 3.1 below, we provide a more detailed analysis of the statistics of the eigenstate coefficients, showing that this is only true for positive integer values of $q \ge 1$. A more precise expression for the full $q$-dependence of the fractal dimensions is given in equation (31) below. Hence, similarly to the GRP model, the intermediate phase of the WRP ensemble is fractal but not multifractal [31, 45, 50] (in which case $D_q$ would actually vary with $q$). As discussed above in section 1, the emergence of such fractal phase is particularly relevant in many contexts.

In summary, the phase diagram of the WRP ensemble is schematically shown in figure 1 and contains three phases: fully delocalized for $\gamma < 1$, Anderson localized for $\gamma > 3/2$, and fractal for $1 < \gamma < 3/2$. The fractal dimensions (for $q \ge 1$) are given by

$$ D_q = D = \begin{cases} 1 & \text{for } \gamma < 1, \\ 3 - 2\gamma & \text{for } 1 < \gamma < 3/2, \\ 0 & \text{for } \gamma > 3/2. \end{cases} \qquad (22) $$

### 3.1. Spectrum of fractal dimensions

In this section we derive the full multifractal spectrum of eigenstate amplitudes using standard perturbation theory. To this end, we closely follow the perturbative calculation done in [31] for the GRP model by adapting it to the WRP ensemble. We denote by
$w_{ij} = |\psi_i(j)|^2$ the squared amplitude of the $i$th eigenvector on site $j$, with $\psi_i(j) = \langle j|\psi_i \rangle$. The first-order correction to the eigenvectors of $A$ reads

$$|\psi_i\rangle = |i\rangle + \sum_{j(\neq i)} \frac{H_{ij}}{a_i - a_j} |j\rangle. \qquad (23)$$

Consequently, for $j \neq i$, the amplitude is

$$w_{ij} = \frac{H_{ij}^2}{(a_i - a_j)^2}. \qquad (24)$$

The convergence of the perturbative series is ensured for $\gamma > 3/2$, where eigenstates are fully localized. For $1 < \gamma < 3/2$, convergence still holds due to the random signs of both $H_{ij}$ and $\Delta_{ij} = a_i - a_j$, as in the RP ensemble.

Equation (24) shows that $w_{ij}$ results from the product of two independent random factors: $x_{ij} = H_{ij}^2$ and $y_{ij} = \Delta_{ij}^{-2}$. The first is the square of a Gaussian variable with variance $\langle H_{ij}^2 \rangle = \nu^2 N^{1-2\gamma}$, while the second has a heavy-tailed distribution $P(y_{ij}) \sim y_{ij}^{-3/2}$ and a typical value of $\mathcal{O}(1)$. Hence, $w_{ij}$ inherits a power-law tail with exponent 3/2 and a typical scale $w_{\text{typ}} \sim N^{1-2\gamma}$. We can thus represent its probability density as

$$P(w_{ij}) = \frac{1}{w_{\text{typ}}} P_{\text{reg}}\left(\frac{w_{ij}}{w_{\text{typ}}}\right) + C \frac{\Theta(w_{ij} > N^{1-2\gamma})}{N^{\gamma-1/2} w_{ij}^{3/2}}, \qquad (25)$$

where $P_{\text{reg}}$ denotes the regular part, and $C$ is fixed by normalization. Imposing the normalization condition $\sum_j |\psi_i(j)|^2 = 1 \Leftrightarrow \langle w_{ij} \rangle = N^{-1}$ introduces an upper cutoff $w_{\text{max}}$ on the singular tail, obtained from

$$N \langle w_{ij} \rangle \sim N^{2-2\gamma} + \frac{C}{N^{\gamma-3/2}} \int_{N^{1-2\gamma}}^{w_{\text{max}}} \mathrm{d}w_{ij} w_{ij}^{-1/2} = 1, \qquad (26)$$

which yields $w_{\text{max}} \sim N^{2(\gamma-3/2)}$. Because amplitudes cannot exceed unity, this expression holds only for $\gamma < 3/2$; when $\gamma > 3/2$, in the localized phase, one must set $w_{\text{max}} = 1$. To restore proper normalization in the localized regime, an additional delta peak at $w_{ij} = 1$ must be included:

$$\hat{P}(w_{ij}) = P(w_{ij}) + A \delta(w_{ij} - 1), \quad \text{for } \gamma > 3/2, \qquad (27)$$

with $A = N^{-1}$, yielding the dominant contribution from the fully localized site.

Next, we define the function $f(\alpha)$ (known as the spectrum of fractal dimensions) in such a way that $\sim N^{f(\alpha)}$ nodes have wave-function amplitudes $|\psi_i(j)|^2 \sim N^{-\alpha}$. For $1 < \gamma < 3/2$, the integration of the tail of $P(w_{ij})$ gives

$$N^{f(\alpha)} = \frac{C}{N^{\gamma-3/2}} \int_{N^{-\alpha}}^{N^{2(\gamma-3/2)}} \mathrm{d}w_{ij} w_{ij}^{-3/2} \sim N^{\alpha/2+3/2-\gamma}, \qquad (28)$$
<figure><img src="image_8.png" /><figcaption>Figure 5. (a) Multifractal spectrum $f(\alpha)$, equation (29), and (b) fractal dimensions $D_q$, equations (31) and (32), for several values of $\gamma$ spanning the three regimes.</figcaption></figure>

valid for $\alpha$ between $\alpha_{min} = 2(3/2 - \gamma)$ and $\alpha_{max} = 2\gamma - 1$. Thus,

$$f(\alpha) = \frac{\alpha}{2} + \frac{3}{2} - \gamma, \quad (\alpha_{\min} < \alpha < \alpha_{\max}). \qquad (29)$$

For $\gamma > 3/2$, the minimum exponent becomes $\alpha_{min} = 0$, and the triangular shape of $f(\alpha)$ at the transition, $f(\alpha) = \alpha/2$, coincides with that of the Anderson model on the Bethe lattice [86, 87].

Alternatively, one may directly compute the moments of the wave-function intensities, $N\langle|\psi_i(j)|^{2q}\rangle \propto N^{-\tau_q}$, using

$$\langle w_{ij}^q \rangle \sim N^{q(1-2\gamma)} + \frac{C}{N^{\gamma-1/2}} \int_{N^{1-2\gamma}}^{N^{2(\gamma-3/2)}} dw_{ij} w_{ij}^{q-3/2}. \qquad (30)$$

For $q < 1/2$, the integral is dominated by typical amplitudes, leading to $\tau_q = q(2\gamma - 1) - 1$, while for $q > 1/2$ large fluctuations dominate, yielding $\tau_q = (q-1)(3-2\gamma)$. The corresponding generalized fractal dimensions, $D_q = \tau_q/(q-1)$, are then ($1 < \gamma < 3/2$)

$$D_q = \begin{cases} 3 - 2\gamma, & q > 1/2, \\ \frac{1 - q(2\gamma - 1)}{1 - q}, & q < 1/2. \end{cases} \qquad (31)$$

In the localized regime ($\gamma > 3/2$), a similar computation gives

$$D_q = \begin{cases} 0, & q > 1/(2\gamma - 1), \\ \frac{1 - q(2\gamma - 1)}{1 - q}, & q < 1/(2\gamma - 1). \end{cases} \qquad (32)$$

The full set of $f(\alpha)$ and $D_q$ curves is displayed in figure 5.

For $\gamma > 3/2$, the distribution $\hat{P}(w_{ij})$ in equation (27) produces a delta peak at $w_{ij} = 1$, resulting in a singular contribution $N^{-1}\delta(w_{ij} - 1)$. The corresponding multifractal spectrum exhibits a spike at $\alpha = 0$; see figure 5(a). In this case, equation 29 must be
modified: one can check that the scaling exponents $\tau_q$ extracted from equation (32) match those of the convex envelope of $f(\alpha)$, namely $f(\alpha) = \alpha/(2\gamma - 1)$ for $0 \le \alpha \le 2\gamma - 1$. An analogous triangular shape with slope smaller than $1/2$ is also observed in the localized regime on random regular graphs [86].

## 4. The cavity method

In this section, we apply the cavity method to determine the diagonal elements of the resolvent matrix of our model (2). This approach not only enables us to re-derive the phase diagram shown in figure 1 beyond perturbation theory, but also provides access to more complex spectral observables, such as the two-point correlations of the energy levels and the spectral compressibility.

The basic idea of this approach is to obtain a self-consistent relation for the resolvent matrix $G = (\lambda_\epsilon \mathbf{1} - \mathbf{H})^{-1}$ of $\mathbf{H}$, which becomes asymptotically exact in the large-$N$ limit. Here, $\lambda_\epsilon = \lambda - i\epsilon$, $\lambda$ being the real energy at which we probe the spectral properties of $\mathbf{H}$, and $\epsilon$ being an imaginary regulator that will be sent to zero at the end of the calculation (after taking the $N \to \infty$ limit). To set the stage, let us assume that we know the resolvent of an $(N-1) \times (N-1)$ matrix $\mathbf{H}^{(1)}$, and let us add an extra row and an extra column (say row 1 and column 1), thereby obtaining an $N \times N$ matrix that we denote $\mathbf{H}$. Using the standard formula of matrix inversion one immediately deduces

$$[G_{11}(\lambda_{\epsilon})]^{-1} = \frac{\operatorname{minor}(\lambda_{\epsilon}\mathbf{1} - \mathbf{H})_{11}}{\det(\lambda_{\epsilon}\mathbf{1} - \mathbf{H})}. \qquad (33)$$

We now use the Schur complement formula (also known as the block matrix inversion formula) to expand the determinant in the denominator in terms of the minor along the first row. After simple algebra one immediately obtains

$$[G_{11}(\lambda_{\epsilon})]^{-1} = \lambda_{\epsilon} - H_{11} - \sum_{i,j=2}^{N} H_{1i} G_{ij}^{(1)}(\lambda_{\epsilon}) H_{j1}, \qquad (34)$$

where the superscript of $G_{ij}^{(1)}$ indicates the element of the resolvent of the $(N-1) \times (N-1)$ matrix $\mathbf{H}^{(1)}$ in the absence of the first row and column, with indices $i, j$ going from 2 to $N$.

This relation is general and exact, without any assumption on the elements $H_{ij}$. However, in order to close these equations and obtain useful relations for the diagonal terms, one must introduce certain approximations. The standard approach consists in neglecting the contributions of the off-diagonal terms $\sum_{i \ne j=2}^{N} H_{1i} G_{ij}^{(1)}(\lambda_{\epsilon}) H_{j1}$, which leads to a closed set of equations for the diagonal elements. However, applying this approximation to the Wishart case (even to the pure Wishart ensemble, without the matrix $\mathbf{A}$) yields an incorrect self-consistent relation—that even fails to reproduce the correct Marčenko–Pastur distribution for the average DoS. As discussed in [2], a different treatment is required for the Wishart ensemble. By substituting $H_{1i}$ and $H_{j1}$ with their
explicit definitions from equations (2) and (5), and separating the resulting double sum over $ℓ$ and $ℓ'$ into contributions with equal and distinct indices, we obtain

$$ \frac{\nu}{M^{\gamma}} \sum_{\ell, \ell'=1}^{M} W_{1\ell} W_{i\ell} W_{1\ell'} W_{j\ell'} = \frac{\nu}{M^{\gamma}} \sum_{\ell=1}^{M} W_{1\ell}^{2} W_{i\ell} W_{j\ell} + \frac{\nu}{M^{\gamma}} \sum_{\ell \neq \ell'} W_{1\ell} W_{i\ell} W_{1\ell'} W_{j\ell'}. \quad (35) $$

In this expression, $W_{1\ell}^2$ is a random variable with mean 1 and variance 2, while $W_{1\ell}W_{1\ell'}$ are Gaussian random variables with zero mean and unit variance. As discussed in [2], for a fixed realization of the elements of $\mathbf{H}^{(1)}$, the left-hand side of the equation above converges in the thermodynamic limit to its average, $\nu M^{-\gamma} \sum_{\ell=1}^{M} W_{i\ell}W_{j\ell} = \nu B_{ij}$. We thus obtain, to the leading order,

$$ \frac{\nu}{M^{\gamma}} \sum_{i,j=2}^{N} G_{ij}^{(1)} \frac{\nu}{M^{\gamma}} \sum_{\ell, \ell'=1}^{M} W_{1\ell} W_{i\ell} W_{1\ell'} W_{j\ell'} \approx \frac{\nu}{M^{\gamma}} \sum_{i,j=2}^{N} G_{ij}^{(1)} (\nu B_{ij} + O(M^{1-\gamma})) \\ = \frac{\nu}{M^{\gamma}} \mathrm{Tr} [\mathbf{G}^{(1)} \nu \mathbf{B}^{(1)}] + O(M^{2(1-\gamma)}). \quad (36) $$

The trace in the expression above can be rewritten as (we omit the superscript (1) to simplify the notation)

$$ \mathrm{Tr} \left[ (\lambda_{\epsilon} \mathbf{1} - \mathbf{A} - \nu \mathbf{B})^{-1} (\nu \mathbf{B} + \mathbf{A} - \lambda_{\epsilon} \mathbf{1} - \mathbf{A} + \lambda_{\epsilon} \mathbf{1}) \right] = \mathrm{Tr} \left[ -\mathbf{1} + (\lambda_{\epsilon} \mathbf{1} - \mathbf{A}) \mathbf{G} \right] \\ = -(N-1) + \sum_{i=2}^{N} (\lambda_{\epsilon} - a_{i}) G_{ii} (\lambda_{\epsilon}). \quad (37) $$

In the large-$N$ limit, the resolvent of the $N \times N$ matrix converges to the resolvent of the $(N-1) \times (N-1)$ matrix (up to corrections of $O(1/N)$). Furthermore, $H_{11} = a_1 + \nu M^{-\gamma} \sum_l W_{1\ell}^2 = a_1 + \nu M^{1-\gamma} + O(M^{1/2-\gamma})$. Neglecting all terms smaller than $M^{1-\gamma}$, one finally obtains the cavity equations for the diagonal elements of the resolvent matrix:

$$ G_{11}^{-1}(\lambda_{\epsilon}) = \lambda_{\epsilon} - a_1 + \nu M^{1-\gamma} \left[ c - 1 - \frac{c}{N} \sum_{i} (\lambda_{\epsilon} - a_i) G_{ii}(\lambda_{\epsilon}) \right]. \quad (38) $$

This equation is asymptotically exact for the WRP ensemble in the large-$N$ limit. In the following, we will use it to extract information on the DoS.

### 4.1. Density of states, local density of states and phase diagram

First, knowledge of the diagonal elements of the resolvent matrix immediately yields the spectral density (8)

$$ \rho^{(e)}(\lambda) = \frac{1}{N\pi} \lim_{\epsilon \to 0} \mathrm{Im} [\mathrm{Tr} \, \mathbf{G}(\lambda_{\epsilon})]. \quad (39) $$
For $\gamma < 1$, the eigenvalues of $\mathbf{H}$ scale as $M^{1-\gamma} \gg 1$, and the support of the DoS grows with $N$. Indeed, upon neglecting the subleading $a_i$ terms, defining $\lambda = M^{1-\gamma} \tilde{\lambda}_\epsilon$ with $\mathrm{Re}(\tilde{\lambda}_\epsilon) \sim \mathcal{O}(1)$, and introducing $g(\lambda_\epsilon) = 1/N \sum_i G_{ii}(\lambda_\epsilon)$, in the large-$N$ limit one obtains [2]

$$ \frac{1}{g(\tilde{\lambda}_{\epsilon})} = \tilde{\lambda}_{\epsilon} + \nu(c-1) - \nu c \tilde{\lambda}_{\epsilon} g(\tilde{\lambda}_{\epsilon}), \qquad (40) $$

from which the Marčenko–Pastur distribution is immediately recovered from equation (39).

In the regime $\gamma > 1$, we have that $M^{1-\gamma} \ll 1$. In this case, it is thus convenient to introduce the small parameter

$$ \eta \equiv \nu M^{1-\gamma} \ll 1, \qquad (41) $$

and the (rescaled) self-energies as

$$ G_{ii}(\lambda_{\epsilon}) = \frac{1}{\lambda_{\epsilon} - a_{i} + \eta \Sigma_{i}(\lambda_{\epsilon})}, \qquad (42) $$

where the $\Sigma_i$'s are of $\mathcal{O}(1)$ (and should not be confused with a summation symbol over $i$). The imaginary part of $G_{ii}$ gives the local density of states (LDoS), representing the contribution of site $i$ to the total DoS:

$$ \rho_i(\lambda) = \sum_{\alpha=1}^{N} |\psi_{\alpha}(i)|^2 \delta(\lambda - \lambda_{\alpha}) = \frac{1}{\pi} \lim_{\epsilon \to 0} \mathrm{Im} G_{ii}(\lambda_{\epsilon}), \qquad (43) $$

so that $\rho^{(e)}(\lambda) = \frac{1}{N} \sum_{i=1}^{N} \rho_i(\lambda)$ (compare with equation (8)). Physically, $\rho_i(\lambda)$ corresponds to the inverse lifetime of a particle created at site $i$ with energy $\lambda$, and provides the order parameter distribution function for Anderson localization.

As shown self-consistently below, the self-energy $\Sigma_i$ turns out to be independent of the index $i$ at leading order for $\gamma > 1$. Plugging the expression (42) into the cavity equation (38) and expanding the terms proportional to $\eta$, one deduces the following self-consistent equation for the self-energy:

$$ \Sigma(\lambda_{\epsilon}) = \frac{\eta c \Sigma(\lambda_{\epsilon})}{N} \sum_{i} \frac{1}{\lambda_{\epsilon} - a_{i}} - 1 + \mathcal{O}(M^{2(1-\gamma)}), \qquad (44) $$

which immediately yields

$$ \Sigma(\lambda_{\epsilon}) \approx - \left( 1 - \frac{\eta c}{N} \sum_{i=1}^{N} \frac{1}{\lambda_{\epsilon} - a_{i}} \right)^{-1}. \qquad (45) $$
In the large-N limit the average of the self-energy over the probability distribution of the diagonal elements can be conveniently rewritten in a more compact form using the definition of the Stieltjes transform of $p_a$:

$$ \mathcal{G}_a(\lambda_\epsilon) = \int da \frac{p_a(a)}{\lambda_\epsilon - a}. \qquad (46) $$

From equation (45) and the definition above, for large N one finds

$$ \Sigma(\lambda_{\epsilon}) \approx -\frac{1}{1 - \eta c \mathcal{G}_a(\lambda_{\epsilon})}, \qquad (47) $$

and from equation (42) one obtains

$$ g(\lambda) = \frac{1}{N} \sum_{i=1}^{N} G_{ii}(\lambda) = \mathcal{G}_a(\lambda + \eta \Sigma) = \mathcal{G}_a \left( \lambda - \frac{\eta}{1 - \eta c \mathcal{G}_a(\lambda)} \right). \qquad (48) $$

Finally, the average DoS in the large-N limit can be written as

$$ \rho(\lambda) = \frac{1}{\pi} \mathrm{Im} \mathcal{G}_a \left( \lambda - \frac{\eta}{1 - \eta c \mathcal{G}_a(\lambda)} \right). \qquad (49) $$

Hence, for $\gamma > 1$, the average DoS of the WRP ensemble coincides with the probability distribution of the diagonal entries, $\rho(\lambda) = p_a(\lambda)$, up to subleading corrections of $\mathcal{O}(\eta)$. In fact, equation (49) follows directly from the Zee formula [83], which expresses the resolvent (i.e. the Stieltjes transform) of the sum of two independent random matrices in terms of their individual resolvents. The same result will be re-derived in section 5.1 using the replica method. We will also return to the finite-N corrections to the average DoS of the WRP ensemble in the different regimes.

The results show that the cavity approach captures the transition of the average DoS taking place at $\gamma = 1$: for $\gamma < 1$, **H** is dominated by **B**, and the average DoS follows, to leading order, the spectral density of **B** (i.e. the Marčenko–Pastur law), up to subleading corrections. Conversely, for $\gamma > 1$, **A** dominates, and the leading contribution to the average DoS is given by the distribution of the diagonal entries, again up to subleading corrections.

We now show that the cavity approach, and in particular equations (42) and (45), also allow one to obtain the second phase transition, occurring at $\gamma = 3/2$, between the intermediate fractal phase and the Anderson-localized one. To this end, we return to equation (45) and introduce the real and imaginary parts of the self-energy as $\Sigma = \sigma - i\hat{\sigma}$ (recall that $\lambda_\epsilon = \lambda - i\epsilon$). For $\gamma > 1$, expanding to leading order in $\eta \ll 1$, one obtains

$$ \begin{aligned} \sigma(\lambda) &\approx -1 - \frac{\eta c}{N} \sum_{i=1}^{N} \frac{1}{\lambda - a_i} + o(\eta), \\ \hat{\sigma}(\lambda) &\approx \frac{\eta c}{N} \sum_{i=1}^{N} \frac{\epsilon}{(\lambda - a_i)^2 + \epsilon^2} + o(\eta). \end{aligned} \qquad (50) $$
Note that we have neglected the imaginary regulator $\epsilon$ with respect to the imaginary part of the self-energy, since the natural scale of the former is of $O(N^{-1})$, while the latter is of order $\eta \propto N^{1-\gamma} \gg N^{-1}$. As discussed in section 3, localization occurs when the corrections to the Green’s function due to the imaginary part of the self-energy (and hence the Thouless energy) are smaller than the mean level spacing. From equation (50) one finds

$$ \hat{\sigma}(\lambda) \approx \eta c \pi p(\lambda). \qquad (51) $$

From equation (42), one immediately obtains the corrections to the imaginary part of the Green’s function due to the perturbation, yielding the Thouless energy $E_T = \eta \hat{\sigma}(\lambda) = \nu^2 c \pi p_a(\lambda) M^{2-2\gamma}$. Anderson localization is then realized when this Thouless energy is smaller than the average gap $1/[Np_a(\lambda)]$, i.e. when

$$ E_T = \nu^2 c \pi M^{2-2\gamma} p_a(\lambda) \ll \frac{1}{N p_a(\lambda)}, \qquad (52) $$

which reproduces the condition $\gamma > 3/2$ obtained from the Fermi golden rule (equation (18)) and Mott’s criterion (equation (16)). The transition to the fully delocalized regime occurs when the Thouless energy exceeds the total bandwidth, i.e. when $\gamma < 1$.

In conclusion, the cavity method allows one to re-derive the phase diagram discussed in section 3 and schematically presented in figure 1, in a more rigorous and controlled way.

## 4.2. Density–density correlations and level compressibility in the intermediate regime

We now apply the cavity approach to compute the two-point density–density correlation function in the intermediate fractal regime, $1 < \gamma < 3/2$, which constitutes one of the main focuses of this work. The two-point function is defined as $\langle \rho(\omega_1) \rho(\omega_2) \rangle_c$. From equations (42), (43) and (50), we first obtain [68]

$$ \langle \rho^{(e)}(\omega_1) \rho^{(e)}(\omega_2) \rangle = \frac{1}{\pi^2 N^2} \sum_{i,j=1}^{N} \langle \rho_i(\omega_1) \rho_j(\omega_2) \rangle, \qquad (53) \\ \rho_i(\omega) = \frac{\eta \hat{\sigma}(\omega)}{(\omega - a_i + \eta \sigma(\omega))^2 + (\eta \hat{\sigma}(\omega))^2}, $$

where $\sigma(\omega)$ and $\hat{\sigma}(\omega)$ are given in equation (50). It is convenient to introduce the discrete version of the Stieltjes transform (46) of $p_a$ at finite $N$ as

$$ S(\lambda_{\epsilon}) = \frac{1}{N} \sum_{i} \frac{1}{\lambda_{\epsilon} - a_{i}} = s(\lambda_{\epsilon}) + i\hat{s}(\lambda_{\epsilon}), \qquad (54) \\ \sigma(\omega) = -1 - \eta c s(\omega), \quad \hat{\sigma}(\omega) = \eta c \hat{s}(\omega). $$

In the intermediate phase, for energy separations of the order of the mean level spacing, $\omega_2 - \omega_1 \sim \Delta \propto N^{-1}$, the level statistics is governed by the universal GOE ensemble. (This
regime actually lies outside the range of validity of the approximations underlying the cavity approach.) For energy separations of the order $\omega_2 - \omega_1 \sim \mathcal{O}(1)$, i.e. much larger than the width of the mini-bands in the local spectrum, the levels are instead expected to be uncorrelated.

The physically relevant regime in the fractal phase is the crossover between RMT behavior at small energy separation and Poisson statistics at large separations, which is expected to occur on the scale of the Thouless energy. For this reason, we consider energies $\omega_1$ and $\omega_2$ separated by intervals of the order of the Thouless energy. For convenience (although this is not strictly necessary), we shift the interval by $\eta$, which slightly simplifies the formulas below. We thus set

$$ \omega_1 = \eta - x\eta^2, \quad \omega_2 = \eta + x\eta^2, \quad x = \mathcal{O}(1). \qquad (55) $$

On these energy scales, one can approximate $s(\omega_j)$ and $\hat{s}(\omega_j)$ (with $j=1$ or 2) as independent of $x$ at the leading order; namely,

$$ s(\omega_j) \approx s(\eta + o(\eta)) \approx s(\eta), \quad \hat{s}(\omega_j) \approx \hat{s}(\eta + o(\eta)) \approx \hat{s}(\eta) \approx \pi [p_a(0) + \eta p'(0)]. \quad (56) $$

Within this approximation, the $\rho_i$'s become random variables that depend only on the set of all $\{a_i\}$'s, and are therefore uncorrelated. One thus obtains

$$ \langle \rho^{(e)}(\omega_1) \rho^{(e)}(\omega_2) \rangle_c \approx \frac{1}{\pi^2 N^2} \sum_{i=1}^{N} \langle \rho_i(\omega_1) \rho_i(\omega_2) \rangle_c = \frac{1}{\pi^2 N} \langle \rho_i(\omega_1) \rho_i(\omega_2) \rangle_c, \quad (57) $$

where the LDoS $\rho_i(\omega)$ is given by the second line of equation (53), and the remaining average on the r.h.s. of equation (57) is intended over $p_a(a)$. Note that $\langle \rho^{(e)}(\omega_1)\rho^{(e)}(\omega_2) \rangle_c$ actually becomes $i$-independent after averaging over $p_a(a)$. This approximation is only valid on the scale of the Thouless energy, and breaks down at the scale of the average spectral gap $1/N$.

The expression above of the density–density correlator allows one to compute the level compressibility $\chi(E)$ introduced in equation (10) using the following relation between the two-point density–density correlation function and the second cumulant of the number of eigenvalues in a given energy window (see e.g. appendix F in [61]):

$$ \kappa_2(\omega_2 - \omega_1) = N \int_{\omega_1}^{\omega_2} d\tilde{\omega}_1 \int_{\omega_1}^{\omega_2} d\tilde{\omega}_2 \langle \rho(\tilde{\omega}_1) \rho(\tilde{\omega}_2) \rangle_c. \qquad (58) $$

Introducing the energy shift $x$ as in equation (55), using equations (53) and (57), and performing the integral over $\tilde{\omega}_1$ and $\tilde{\omega}_2$ (with the change of variable $\tilde{\omega}_1 = \eta + x\eta^2$ and $\tilde{\omega}_2 = \eta - x\eta^2$) before the one over the $a_i$'s, one finally obtains

$$ \pi^2 \kappa_2(x) = \int da p_a(a) \left[ \begin{aligned} & \operatorname{atan}\left( \frac{\eta^2 x + a + \eta^2 s(\eta)}{\eta^2 c \hat{s}(\eta)} \right) - \operatorname{atan}\left( \frac{a - \eta^2 x + \eta^2 s(\eta)}{\eta^2 c \hat{s}(\eta)} \right) \end{aligned} \right]^2 \\ & - \left\{ \int da p_a(a) \left[ \begin{aligned} & \operatorname{atan}\left( \frac{\eta^2 x + a + \eta^2 s(\eta)}{\eta^2 c \hat{s}(\eta)} \right) - \operatorname{atan}\left( \frac{a - \eta^2 x + \eta^2 s(\eta)}{\eta^2 c \hat{s}(\eta)} \right) \end{aligned} \right] \right\}^2. \qquad (59) $$
Since the difference of the two atan's is non-zero only in a small interval centered around zero and of width $\eta^2$, one can perform the change of variable $a = \eta^2 \tilde{a}$ and replace $p_a(\eta^2 \tilde{a})$ by $p_a(0)$ in the integrals above. The disconnected part of $\kappa_2$, corresponding to the integral in the second line, gives $\kappa_1^2(x) = [2\eta^2 p_a(0)(x + s(\eta))]^2$. The term in the second line is thus of $O(\eta^4)$, and can be neglected with respect to the term in the first line, which is of $O(\eta^2)$. For a symmetric distribution $p_a$, the real part of the Stieltjes transform (46) computed at energy $\eta$ is of $O(\eta)$, and can also be neglected at the leading order. Furthermore, it is convenient to change the variable again as $\tilde{a} \to c\hat{s}\tilde{a}$. Since $\hat{s} \approx \pi p_a(0)$, we have

$$ \kappa_2(x) \approx \frac{\eta^2 p_a^2(0) c}{\pi} \int_{-\infty}^{+\infty} d\tilde{a} \left[ \operatorname{atan}\left(\frac{x}{c\pi p_a(0)} + \tilde{a}\right) - \operatorname{atan}\left(\tilde{a} - \frac{x}{c\pi p_a(0)}\right) \right]^2 . \quad (60) $$

Plugging this result into equation (10) and dividing by the first cumulant $\kappa_1(x) \approx 2\eta^2 p_a(0)x$, one finally obtains the spectral compressibility on the scale of the Thouless energy:

$$ \chi(x) \approx \frac{c p_a(0)}{2\pi x} \int_{-\infty}^{+\infty} d\tilde{a} \left[ \operatorname{atan}\left(\frac{x}{c\pi p_a(0)} + \tilde{a}\right) - \operatorname{atan}\left(\tilde{a} - \frac{x}{c\pi p_a(0)}\right) \right]^2 . \quad (61) $$

Using the property (see e.g. 1.625 in [88])

$$ \operatorname{atan}(a) - \operatorname{atan}(b) = \operatorname{atan}\left(\frac{a-b}{1+ab}\right), \quad (62) $$

upon calling $y \equiv \frac{x}{c\pi p_a(0)}$ in equation (61) and changing variables in the integral as $u = \tilde{a}\sqrt{1+y^2}$, we obtain

$$ \chi \left( y = \frac{\omega_2 - \omega_1}{2\pi c p_a(0) \eta} \right) = \frac{\sqrt{1+y^2}}{\pi^2 y} \int_0^\infty du \left\{ \operatorname{atan} \left[ \frac{2y}{u^2(1+y^2) + 1 - y^2} \right] \right\}^2, \quad (63) $$

where we stress that we have chosen the branch atan($z$) $\in [0, \pi]$. Upon integrating by parts and performing some algebra [88, 89], the integral over $u$ can be computed explicitly to give equation (13), which is exactly the same universal scaling function that some of us found for the GRP ensemble (see equations (126) and (127) of [61]). Alternatively (and equivalently), in appendix A we provide a derivation of the scaling function (13) starting from equation (61). In particular, the asymptotics of the function $\chi(y)$ can be checked to give

$$ \chi(y) \approx \begin{cases} \frac{y}{\pi}, & y \ll 1, \\ 1 - \frac{2(1 + \ln y)}{\pi y}, & y \gg 1, \end{cases} \quad (64) $$

showing that $\chi(y)$ interpolates between Wigner-Dyson statistics at low energy, and Poisson statistics at higher energy. This function is plotted as a dot-dashed line in figures 3 and 4.
# 5. The replica method

In this section we derive the average DoS reported in equation (49) using the replica method [90]. We also apply the replica strategy to calculate the level compressibility given in equation (13). In passing, we derive in section 5.2 the full counting statistics of the WRP model.

## 5.1. Density of states

We report here the main steps of the derivation, while we defer its details to appendix B. Using the Edwards–Jones formula [91], we first express the average DoS as

$$ \rho(\lambda) = -\frac{2}{N\pi} \lim_{\epsilon \to 0^+} \mathrm{Im} \frac{\partial}{\partial \lambda} \langle \ln \mathcal{Z}(\lambda - i\epsilon) \rangle, \qquad (65) $$

where we introduced the partition function

$$ \mathcal{Z}(\lambda) = \int_{\mathbb{R}^N} d^N r \, e^{-\frac{1}{2} r^T (\lambda 1 - \mathbf{H}) r}. \qquad (66) $$

The average of the logarithm in equation (65) can be expressed using the replica trick as

$$ \langle \ln \mathcal{Z}(\lambda) \rangle = \lim_{n \to 0} \frac{1}{n} \ln \langle \mathcal{Z}^n(\lambda) \rangle. \qquad (67) $$

Using standard techniques, one can then recast the average as

$$ \langle \mathcal{Z}^n(\lambda) \rangle = \int \mathcal{D}\phi \mathcal{D}\hat{\phi} \mathcal{D}\psi \mathcal{D}\hat{\psi} \exp \left\{ -\sqrt{NM} \mathcal{S}_n \left[ \phi, \hat{\phi}, \psi, \hat{\psi}; \lambda \right] \right\}, \qquad (68) $$

with the action $\mathcal{S}_n[\phi, \hat{\phi}, \psi, \hat{\psi}; \lambda]$ given by

$$ \begin{aligned} \mathcal{S}_n[\phi, \hat{\phi}, \psi, \hat{\psi}; \lambda] = & \frac{i}{\sqrt{c}} \int d\vec{u} \, \hat{\phi}(\vec{u}) \phi(\vec{u}) + i\sqrt{c} \int d\vec{r} \, \hat{\psi}(\vec{r}) \psi(\vec{r}) \\ & - \sqrt{c} \ln \int d\vec{r} \, \varphi_a(-\vec{r}^2/2) \exp \left[ -\frac{i}{2} \lambda \vec{r}^2 + i \hat{\psi}(\vec{r}) \right] \\ & - \frac{1}{\sqrt{c}} \ln \int d\vec{u} \, \exp \left[ -\frac{1}{2} \vec{u}^2 + i \hat{\phi}(\vec{u}) \right] + \frac{i}{2} \sqrt{c\eta} \int d\vec{r} d\vec{u} \, \phi(\vec{u}) \psi(\vec{r}) (\vec{u} \cdot \vec{r})^2, \end{aligned} \qquad (69) $$

where $\eta = \nu M^{(1-\gamma)}$ and $\varphi_a$ is the characteristic function of $p_a$, i.e.

$$ \varphi_a(y) = \int da \, p_a(a) e^{-iay}. \qquad (70) $$

(Here and in appendix B, we denote by $\mathbf{r}$ a vector in $\mathbb{R}^N$, and by $\vec{r}$ a vector in the replica space $\mathbb{R}^n$.) Next, using saddle-point evaluation, we get that $\langle \mathcal{Z}^n(\lambda) \rangle \approx$
$\exp(-\sqrt{NM}\mathcal{S}_n^{\text{sp}})$, where $\mathcal{S}_n^{\text{sp}}$ denotes the action evaluated at the saddle point. Combining this with equations (65) and (67), we finally obtain the DoS as

$$ \rho(\lambda) = \frac{1}{\pi} \lim_{\epsilon \to 0^+} \text{Im} G(\lambda_{\epsilon}) + \mathcal{O}(1/N), \qquad (71) $$

where $G(\lambda) \propto \lim_{n \to 0} \partial_{\lambda} \mathcal{S}_n^{\text{sp}}$ (which actually coincides with the resolvent associated with the WRP model; compare with equation (39)) satisfies the self-consistent equation

$$ G(\lambda) = i \int_{0}^{\infty} dz \, \varphi_a(-z) \exp \left[ -iz \left( \lambda - \frac{\eta}{1 - c\eta G(\lambda)} \right) \right]. \qquad (72) $$

Equivalently, upon rewriting the characteristic function as in equation (70), switching the order of the two integrals, computing the integral over $z$ first and then the one over $a$, we obtain

$$ G(\lambda) = \mathcal{G}_a \left( \lambda - \frac{\eta}{1 - c\eta G(\lambda)} \right), \qquad (73) $$

where $\mathcal{G}_a$ is the Stieltjes transform of $p_a$, defined as in equation (46). This is consistent with the result found with the cavity method in equation (49), since $G(\lambda) = \mathcal{G}_a(\lambda) + \mathcal{O}(\eta)$. We also note that this result is equivalent to the Zee formula [83], which allows one to compute the resolvent $G_{1+2}$ of the sum of two mutually free random matrices, given their respective resolvents $G_1$ and $G_2$. Indeed, the Zee formula can be written as a self-consistent equation (see e.g. appendix C in [61])

$$ G_{1+2}(z) = G_1(z - R_2(G_{1+2}(z))), \qquad (74) $$

in terms of the $R$-transform of the second matrix, i.e. $R_2$. In our case, the $R$-transform of Wishart matrices is $R_W(z) = 1/(1 - cz)$ [2]. Now, using the fact that rescaling a random matrix by a factor $\eta$ scales its $R$-transform as $R_{\eta W}(z) = \eta R_W(\eta z)$, and plugging this into equation (74), we again obtain equation (73).

In the case of Cauchy-distributed $p_a$ (centered at $\mu$ and of width $\omega$), the Stieltjes transform is known in closed form: $\mathcal{G}_{\text{Cauchy}}(\lambda) = 1/(\lambda - \mu \pm i\omega)$, where the $\pm$ branches correspond to $\text{Im}\lambda > 0$ or $\text{Im}\lambda < 0$, respectively. This renders equation (73) a quadratic equation that can be easily solved (however, note that any other choice of $p_a$ could be evaluated numerically using equation (72)). The average spectral density in the particular case of Cauchy-distributed $p_a$ is presented in figure 6, where it is compared to the numerical diagonalization of samples of the WRP ensemble.

## 5.2. Full counting statistics

Here, we apply the replica method to derive the level compressibility given in equation (13). To this end, we follow the procedure introduced in [92–94], and more recently applied in [61], to the case of the GRP model. Accordingly, we first obtain the cumulant generating function of the number of eigenvalues in a given interval (see equation (9)) in the large-$N$ limit—this quantity is also known as full counting statistics.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** DoS of the WRP ensemble for a Cauchy-distributed $p_a$ in (a) the intermediate phase ($1 < \gamma = 5/4 < 3/2$), and in (b) the delocalized phase ($\gamma = 3/4 < 1$). The solid lines correspond to our prediction (73), while the histograms are obtained from exact diagonalization of $10^5$ samples of the WRP ensemble, with $c = 1/2$, $\nu = 1$, and varying the matrix size $N$. Upon increasing $N$, the densities converge either to $p_a$ (which is the Cauchy distribution here) in (a), or to the Marčenko-Pastur distribution in (b)—this corresponds to the transition sketched in figure 1. Note that, in the delocalized phase, we have rescaled the eigenvalues $\lambda$ by $\eta$ in order to compensate for the growing support.</figcaption></figure>

From this expression, the level compressibility can then be retrieved using equation (10). Here, we will mainly focus on the intermediate phase ($1 < \gamma < 3/2$), and on intervals of the scale of the Thouless energy.

The number of eigenvalues in an interval delimited by $\alpha$ and $\beta$ is given by

$$I_N[\alpha, \beta] = \sum_{i=1}^{N} [\Theta(\beta - x_i) - \Theta(\alpha - x_i)], \qquad (75)$$

where $\Theta(x)$ is the Heaviside distribution. Using the standard identity

$$\Theta(-x) = \frac{1}{2\pi i} \lim_{\epsilon \to 0^+} [\ln(x + i\epsilon) - \ln(x - i\epsilon)], \qquad (76)$$

from complex analysis, one can show that

$$I_N[\alpha, \beta] = -\frac{1}{\pi i} \lim_{\epsilon \to 0^+} \ln \left[ \frac{\mathcal{Z}(\beta - i\epsilon) \mathcal{Z}(\alpha + i\epsilon)}{\mathcal{Z}(\beta + i\epsilon) \mathcal{Z}(\alpha - i\epsilon)} \right], \qquad (77)$$

where $\mathcal{Z}$ is the partition function appearing in the calculation of the spectral density in equation (66) (note that, for $\mathcal{Z}(\alpha + i\epsilon)$ and $\mathcal{Z}(\beta + i\epsilon)$, one actually has to compute $\int_{\mathbb{R}^N} d^N r e^{\frac{i}{2} r^T (\lambda \mathbf{1} - \mathbf{H}) r}$ instead of $\int_{\mathbb{R}^N} d^N r e^{-\frac{i}{2} r^T (\lambda \mathbf{1} - \mathbf{H}) r}$ to ensure convergence of the integral).
Now, assuming that one can exchange the limit $\epsilon \to 0^+$ and the logarithm, one writes the cumulant generating function of the random variable $I_N[\alpha, \beta]$ as

$$ \mathcal{F}_{[\alpha, \beta]}(s) \equiv \ln \langle e^{-sI_N[\alpha, \beta]} \rangle = \lim_{\epsilon \to 0^+} \ln \left\langle \left[ \mathcal{Z}(\beta_{\epsilon}^*) \mathcal{Z}(\alpha_{\epsilon}) \right]^{is/\pi} \left[ \mathcal{Z}(\beta_{\epsilon}) \mathcal{Z}(\alpha_{\epsilon}^*) \right]^{-is/\pi} \right\rangle . \quad (78) $$

Once again, in the spirit of the replica method, we first calculate

$$ Q_{[\alpha, \beta]}(n_{\pm}) = \left\langle \left[ \mathcal{Z}(\beta_{\epsilon}^*) \mathcal{Z}(\alpha_{\epsilon}) \right]^{n_{\pm}} \left[ \mathcal{Z}(\beta_{\epsilon}) \mathcal{Z}(\alpha_{\epsilon}^*) \right]^{n_{-}} \right\rangle, \quad (79) $$

where $n_{\pm}$ are two independent integers. Next, we perform the analytic continuation of $n_{\pm}$ to the imaginary axis to express $\mathcal{F}_{[\alpha, \beta]}(s)$ in the form

$$ \mathcal{F}_{[\alpha, \beta]}(s) = \lim_{\epsilon \to 0^+} \ln \lim_{n_{\pm} \to \pm is/\pi} Q_{[\alpha, \beta]}(n_{\pm}). \quad (80) $$

Following steps similar to those used in the calculation of the DoS, we obtain (see appendix B)

$$ Q_{[\alpha, \beta]}(n_{\pm}) \propto \int \mathcal{D}\phi \mathcal{D}\hat{\phi} \mathcal{D}\psi \mathcal{D}\hat{\psi} \exp \left\{ -\sqrt{NM} \mathcal{S}_{n_{\pm}} \left[ \phi, \hat{\phi}, \psi, \hat{\psi}; \hat{\Lambda} \right] \right\}, \quad (81) $$

with the action given by

$$ \mathcal{S}_{n_{\pm}}[\phi, \hat{\phi}, \psi, \hat{\psi}; \hat{\Lambda}] = \frac{i}{\sqrt{c}} \int d\vec{u} \, \hat{\phi}(\vec{u}) \phi(\vec{u}) + i\sqrt{c} \int d\vec{r} \, \hat{\psi}(\vec{r}) \psi(\vec{r}) \\ - \sqrt{c} \ln \int d\vec{r} \, \varphi_a \left(-\frac{1}{2}\vec{r}\hat{L}\vec{r}\right) \exp\left[-\frac{i}{2}\vec{r}\hat{\Lambda}\vec{r} + i\hat{\psi}(\vec{r})\right] \\ - \frac{1}{\sqrt{c}} \ln \int d\vec{u} \exp\left[-\frac{1}{2}\vec{u}^2 + i\hat{\phi}(\vec{u})\right] + \frac{i}{2}\sqrt{c\eta} \int d\vec{r} d\vec{u} \, \phi(\vec{u})\psi(\vec{r}) (\vec{u}\hat{L}\vec{r})^2, \quad (82) $$

and where we introduced the block matrices

$$ \hat{\Lambda} = \begin{pmatrix} \alpha_{\epsilon} \mathbf{1}_{n_{+}} & & \\ & -\beta_{\epsilon}^{*} \mathbf{1}_{n_{+}} & \\ & & \beta_{\epsilon} \mathbf{1}_{n_{-}} \\ & & & -\alpha_{\epsilon}^{*} \mathbf{1}_{n_{-}} \end{pmatrix}, \quad \hat{L} = \begin{pmatrix} \mathbf{1}_{n_{+}} & & \\ & -\mathbf{1}_{n_{+}} & \\ & & \mathbf{1}_{n_{-}} \\ & & & -\mathbf{1}_{n_{-}} \end{pmatrix}. \quad (83) $$

Next, using the saddle-point approximation, one can obtain a re-parametrisation of the action in terms of two matrices $\hat{K}$ and $\hat{C}^{-1}$ (see appendix B):

$$ \mathcal{S}_{n_{\pm}}[\hat{K}, \hat{C}^{-1}; \hat{\Lambda}] = -\sqrt{c} \ln \left[ \int_{-\infty}^{\infty} da \, p_a(a) \exp \left\{ -\frac{1}{2} \ln \det (\hat{C}^{-1} - ia\hat{L}) \right\} \right] \\ + \frac{1}{2\sqrt{c}} \mathrm{Tr} \log (\hat{1} + iq\eta\hat{K}) - \frac{i}{2}\sqrt{c\eta} \mathrm{Tr} \left[ \hat{K} (\hat{1} + iq\eta\hat{K})^{-1} \right], \quad (84) $$
where $\hat{C}^{-1}$ and $\hat{K}$ have to satisfy the self-consistent relations

$$ \hat{C}^{-1} = i\hat{\Lambda} + i\eta\hat{L} (\hat{1} + iq\eta\hat{K})^{-1}\hat{L}, \quad \text{and} \quad \hat{K} = -i\hat{L}\mathcal{G}_a \left( (i\hat{L}\hat{C})^{-1} \right), \quad (85) $$

where $\mathcal{G}_a$ was given in equation (46). We now look for a block diagonal solution of equation (85), i.e. of the form

$$ \hat{K} \equiv \begin{pmatrix} k_{\alpha} \mathbb{1}_{n_{+}} & \bar{k}_{\beta} \mathbb{1}_{n_{+}} & \mathbb{1}_{n_{-}} & \bar{\mathbb{1}} \\ \mathbb{1}_{n_{+}} & k_{\beta} \mathbb{1}_{n_{-}} & \mathbb{1}_{n_{-}} & \bar{k}_{\alpha} \mathbb{1}_{n_{-}} \\ \mathbb{1}_{n_{+}} & \mathbb{1}_{n_{-}} & \bar{k}_{\alpha} \mathbb{1}_{n_{-}} & \mathbb{1}_{n_{-}} \end{pmatrix}, \quad \hat{C}^{-1} \equiv \begin{pmatrix} \Delta_{\alpha}^{-1} \mathbb{1}_{n_{+}} & \bar{\Delta}_{\beta}^{-1} \mathbb{1}_{n_{+}} & \mathbb{1}_{n_{-}} & \bar{\mathbb{1}} \\ \mathbb{1}_{n_{+}} & \Delta_{\beta}^{-1} \mathbb{1}_{n_{+}} & \mathbb{1}_{n_{-}} & \bar{\Delta}_{\alpha}^{-1} \mathbb{1}_{n_{-}} \end{pmatrix}. \quad (86) $$

This Ansatz can then be inserted into equation (84) to derive

$$ \begin{aligned} S_{n_{\pm}} = & -\sqrt{c} \ln \left\{ \int_{-\infty}^{\infty} da \, p_a(a) e^{-\frac{n_{+}}{2} \ln [(\Delta_{\alpha}^{-1} - ia)(\bar{\Delta}_{\beta}^{-1} + ia)] - \frac{n_{-}}{2} \ln [(\Delta_{\beta}^{-1} - ia)(\bar{\Delta}_{\alpha}^{-1} + ia)]} \right\} \\ & + \frac{1}{2\sqrt{c}} \left\{ n_{+} \ln [(1 + i\eta ck_{\alpha})(1 + i\eta c\bar{k}_{\beta})] + n_{-} \ln [(1 + i\eta ck_{\beta})(1 + i\eta c\bar{k}_{\alpha})] \right\} \\ & - \frac{i}{2} \sqrt{c\eta} \left[ n_{+} \left( \frac{k_{\alpha}}{1 + i\eta ck_{\alpha}} + \frac{\bar{k}_{\beta}}{1 + i\eta c\bar{k}_{\beta}} \right) + n_{-} \left( \frac{k_{\beta}}{1 + i\eta ck_{\beta}} + \frac{\bar{k}_{\alpha}}{1 + i\eta c\bar{k}_{\alpha}} \right) \right]. \end{aligned} \quad (87) $$

Now, taking the limit $n_{\pm} \to \pm is/\pi$, we obtain

$$ \begin{aligned} S_{\pm \frac{is}{\pi}} = & -\sqrt{c} \ln \left\{ \int_{-\infty}^{\infty} da \, p_a(a) \exp \left( -\frac{is}{2\pi} \ln \left[ \frac{(\Delta_{\alpha}^{-1} - ia)(\bar{\Delta}_{\beta}^{-1} + ia)}{(\Delta_{\beta}^{-1} - ia)(\bar{\Delta}_{\alpha}^{-1} + ia)} \right] \right) \right\} \\ & + \frac{is}{2\pi\sqrt{c}} \ln \left[ \frac{(1 + i\eta ck_{\alpha})(1 + i\eta c\bar{k}_{\beta})}{(1 + i\eta ck_{\beta})(1 + i\eta c\bar{k}_{\alpha})} \right] \\ & - \frac{\eta s}{2\pi\sqrt{c}} \left[ \left( \frac{k_{\alpha}}{1 + i\eta ck_{\alpha}} + \frac{\bar{k}_{\beta}}{1 + i\eta c\bar{k}_{\beta}} \right) - \left( \frac{k_{\beta}}{1 + i\eta ck_{\beta}} + \frac{\bar{k}_{\alpha}}{1 + i\eta c\bar{k}_{\alpha}} \right) \right]. \end{aligned} \quad (88) $$

This finally allows us to obtain the cumulant generating function $\mathcal{F}_{[\alpha, \beta]}(s)$, which follows from equations (80) and (81) as

$$ \mathcal{F}_{[\alpha, \beta]}(s) = \sqrt{NM} \lim_{\epsilon \to 0^+} S_{\pm is/\pi} + \mathcal{O}(N^{-\gamma}). \quad (89) $$

This concludes the replica calculation. To get the moments of the number of eigenvalues in an interval $I_N[\alpha, \beta]$, one then proceeds as follows:

(i) The self-consistent equations (85) can be solved numerically to find the eight elements of the block matrices $\hat{K}$ and $\hat{C}^{-1}$ in equation (86). Note that these are actually eight nonlinear equations, which can be solved with minimal numerical effort for any reasonable choice of $p_a$.
(ii) Inserting these matrix elements into equation (88) makes the cumulant generating function fully explicit, in spite of its seemingly complicated integral form.

(iii) Expanding equation (88) in powers of s and using equation (89), one can identify the cumulants $\kappa_j$ as

$$ \kappa_j[\alpha, \beta] = (-1)^j \partial_s^j \mathcal{F}_{[\alpha, \beta]}(s) \Big|_{s=0} \qquad (90) $$

We conclude by pointing out that equation (77), i.e. the starting point of our calculation, was actually obtained by adopting the identity $\ln(ab) = \ln a + \ln b$, which is, however, not satisfied in general by the complex logarithm (whose principal branch is bounded within $(-\pi, \pi]$ [95]). This issue is usually (and quite remarkably) solved via the introduction of replicas, thanks to the so-called *folding-unfolding mechanism* [94]. However, the cumulant generating function obtained in equations (88) and (89) cannot be immediately recognized as a real quantity, as one would have hoped for in general. In other random matrix ensembles, for which the spectral density and/or the chosen interval are symmetric around the origin, the vanishing of the imaginary part of $\mathcal{F}_{[\alpha, \beta]}(s)$ can actually be proven analytically [61, 93]. In the present case, the asymmetry of the spectral density prevented us from carrying out such a proof; however, we have checked that indeed $\mathcal{F}_{[\alpha, \beta]}(s)$ becomes real when $N \to \infty$ for $\gamma > 1$, in which case the DoS also becomes symmetric, and we have considered the real part of $\mathcal{F}_{[\alpha, \beta]}(s)$ otherwise$^8$.

### 5.3. Level compressibility

By specializing the interval to $[\alpha, \beta] = [-E + \eta, E + \eta]$, we can now focus on the ratio between the first two cumulants given by equation (90), which corresponds to the level compressibility $\chi(E)$ introduced in equation (10). The result is plotted in figure 3, where it is tested against numerical diagonalization of large sample random matrices. We can generically distinguish three regimes:

(i) For energies $E < \Delta$, where $\Delta \approx 1/N$ is the mean level spacing, the saddle-point calculation breaks down—as expected, since in our derivation we have treated the eigenvalue density as a continuous distribution. At these energies, the eigenvalue statistics is dominated by level repulsion, as is typically the case in the RMT regime. Here, $\chi(E)$ turns out to be well approximated by the level compressibility of a GOE matrix (see e.g. appendix E.3 in [61]):

$$ \chi_{\text{GOE}}(y) = \frac{1}{2\pi^2 y} \left\{ \begin{aligned} & [\text{Si}(2\pi y)]^2 - 2\text{Ci}(4\pi y) - \pi \text{Si}(2\pi y) \\ & + 2[-4\pi y \text{Si}(4\pi y) + 2\pi^2 y + \log(4\pi y) - \cos(4\pi y) + \gamma_E + 1] \end{aligned} \right\}, \quad (91) $$

<footnote>$^8$ The method adopted here, although with a different Ansatz, was applied in [94] to characterize the *index* (i.e. the number of eigenvalues in the interval $(-\infty, \beta]$) within the diluted Wishart ensemble, whose spectrum is in fact not symmetric. Inspecting this quantity within the WRP ensemble would be insightful in the future, in view of better assessing the limitations of the replica method.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparison between the universal scaling form (13) of the level compressibility $\chi_T$ (solid line), and the analytical replica solution for $\chi(E)$ (symbols, see section 5.2), computed for intervals of the order of the Thouless energy. Upon increasing the matrix size $N$—or equivalently, upon decreasing $\eta \propto N^{1-\gamma}$ for $\gamma > 1$—the predicted level compressibility approaches the universal form $\chi_T$, as conjectured in [61].</figcaption></figure>

where $Ci(z) = -\int_z^\infty \cos(t)/t dt$ and $Si(z) = \int_0^z \sin(t)/t dt$ are the cosine-integral and sine-integral functions, respectively, while $\gamma_E$ is the Euler–Mascheroni constant. The fact that $\chi_{GOE}(y)$ also well describes the level compressibility of a Wishart matrix is best rationalized within the Coulomb gas interpretation [1]: indeed, at these scales, correlations between eigenvalues originate from the Coulomb gas interaction, which is the same for both GOE and real Wishart matrices.

(ii) Around the Thouless energy, i.e. for $E \sim E_T \propto N^{2(1-\gamma)}$, the level compressibility is well described by the universal scaling form (13) found in [61] within the GRP model. The agreement with the numerical results improves upon increasing $N$, as we exemplified in figure 4(a). Moreover, the independence of $\chi(E \sim E_T)$ from the ratio $c = N/M$, which is an expected but nontrivial feature, is demonstrated in figure 4(b). Finally, in figure 7 we show that the replica prediction also approaches $\chi_T$ upon increasing $N$ (i.e. upon decreasing $\eta \propto N^{1-\gamma}$).

(iii) For energies of $O(1)$, the eigenvalues behave as uncorrelated random variables (equation (11)), whence (see e.g. appendix A in [61])

$$ \chi(E) \sim \chi_{\text{iid}}(E) = 1 - \frac{\langle I_N[-E, E] \rangle}{N} = 1 - \int_{-E}^{E} d\lambda p_a(\lambda). \qquad (92) $$

## 6. Dyson Brownian motion

In this section, we derive the phase diagram of the WRP model using yet another approach, similar to the one first introduced in [96, 97], and later applied to the GRP model in [33]. The idea is to interpret equation (2) as a matrix-valued stochastic process,
for which **A** represents the initial condition, and **B** a perturbation that is turned on at the fictitious 'time' $t=0$; the matrix **H** is eventually recovered at the final time $t=T$. Under the effect of this perturbation, the eigenvalues and eigenvectors themselves become stochastic processes, and inspecting their behavior can give us information on the limiting eigenvalue density and the statistics of the phase (localized, delocalized or fractal).

We start by setting up the problem. We study the matrix $\mathbf{H}(t) = \mathbf{A} + \mathbf{X}(t)\mathbf{X}(t)^T$, where $\mathbf{X}(t)$ is an $N \times M$ matrix-valued Brownian motion, and $\mathbf{A}$ is the same diagonal random matrix as in equation (2). At each time step, we can write the evolution of $\mathbf{X}(t)$ as $\mathbf{X}(t+dt) = \mathbf{X}(t) + \mathbf{g}(t)$, with $\mathbf{g}(t)$ being an i.i.d. uncorrelated $N \times M$ Gaussian noise with zero mean and variance proportional to $dt$, i.e. $\langle g_{ik}(t) \rangle = 0$ and $\langle g_{ik}(t)g_{jl}(t') \rangle = \sigma^2 dt \delta_{ij}\delta_{kl} \delta(t-t')$. The stochastic evolution equation of $\mathbf{H}(t)$ then follows as

$$ 
\begin{aligned}
\mathbf{H}(t+dt) &= \mathbf{A} + \mathbf{X}(t+dt)\mathbf{X}(t+dt)^T \\
&= \mathbf{A} + \mathbf{X}(t)\mathbf{X}(t)^T + \mathbf{g}(t)\mathbf{X}(t)^T + \mathbf{X}(t)\mathbf{g}(t)^T + \mathbf{g}(t)\mathbf{g}(t)^T \\
&= \mathbf{H}(t) + \delta\mathbf{H}(t),
\end{aligned} 
\quad (93) $$

where in the last line we defined $\delta\mathbf{H}(t) = \mathbf{g}(t)\mathbf{X}(t)^T + \mathbf{X}(t)\mathbf{g}(t)^T + \mathbf{g}(t)\mathbf{g}(t)^T$. Upon requiring that, at time $T$, $\mathbf{H}(t=T) = \mathbf{A} + \nu M^{-\gamma} \mathbf{W}\mathbf{W}^T$, and by using that $\sum_{l=1}^M X_{il}X_{jl} \sim Mt$, we get that the stopping time of this process must be given by $T = \nu M^{-\gamma}$. Now, applying perturbation theory up to second order, one can derive the following stochastic differential equations for the eigenvalues $\lambda_i(t)$ and the components of the associated eigenvectors $\langle n|\psi_i(t) \rangle = \psi_i(n;t)$ of $\mathbf{H}(t)$, starting from the initial condition $\lambda_i(t=0) = a_i$ and $\psi_i(n;t=0) = \delta_{in}$ (we omitted the time dependencies for clarity, and the details of the derivation are given in appendix C):

$$ \frac{d\lambda_i}{dt} = \sigma^2 \left[ M + \sum_{j \neq i} \frac{\lambda_i + \lambda_j - \langle a \rangle_i - \langle a \rangle_j}{\lambda_i - \lambda_j} \right] + 2\sigma \sqrt{\lambda_i - \langle a \rangle_i} \zeta_i, \quad (94) $$

and

$$ \frac{d\psi_i(n)}{dt} = \sigma^2 \sum_{j \neq i} \psi_j(n) \frac{\langle a \rangle_{ij}}{\lambda_i - \lambda_j} \left( \frac{1}{\lambda_i - \lambda_j} - \sum_{l \neq i} \frac{1}{\lambda_i - \lambda_l} \right) \\ - \frac{1}{2} \sigma^2 \psi_i(n) \sum_{j \neq i} \frac{\lambda_i + \lambda_j - \langle a \rangle_i - \langle a \rangle_j}{\lambda_i - \lambda_j} + \sigma \sum_{j \neq i} \psi_j(n) \frac{\zeta_{ij}}{\lambda_i - \lambda_j}. \quad (95) $$

Here, $\zeta_i$ is a delta-correlated Gaussian white noise, while $\zeta_{ij}$ is a correlated Gaussian noise with zero mean and variance:

$$ \langle \zeta_{ij}(t) \zeta_{kl}(t') \rangle = \delta(t-t') [(\lambda_i + \lambda_j)(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) \\ - (\delta_{ik}\langle a \rangle_{jl} + \delta_{il}\langle a \rangle_{jk} + \delta_{jk}\langle a \rangle_{il} + \delta_{jl}\langle a \rangle_{ik})]. \quad (96) $$

We also defined

$$ \langle a \rangle_i(t) = \langle \psi_i(t) | \mathbf{A} | \psi_i(t) \rangle = \sum_n a_n |\psi_i(n;t)|^2, \quad (97) $$
$$ \langle a \rangle_{ik}(t) = \langle \psi_i(t) | \mathbf{A} | \psi_k(t) \rangle = \sum_n a_n \psi_i(n; t) \psi_k(n; t), \quad (98) $$

which implies that, in this model, eigenvalues are coupled to eigenvectors (this was not the case in the GRP model—see equations (12) and (13) in [33]). Moreover, one can note that setting A to 0 correctly renders the Dyson Brownian motion equation for the eigenvalues and eigenvectors of a pure Wishart process [98]. Although in general these equations are complicated to solve, here we will merely be concerned with their short-time behavior, because the stopping time $T \sim N^{-\gamma}$ is vanishingly small for all values of $\gamma$ of interest.

We begin by considering the evolution of the eigenvalues in equation (94). At short times, only the first term in equation (94) gives a relevant contribution, i.e. $\frac{d\lambda_i}{dt} \simeq \sigma^2 M$, because $\lambda_i(t=0) = a_i = \langle a \rangle_i(t=0)$. Thus, after a time $T \sim N^{-\gamma} \ll 1$, the eigenvalues have moved by

$$ \delta\lambda_i(T) = \lambda_i(t) - a_i \sim \sigma^2 MT \sim N^{1-\gamma}. \quad (99) $$

This means that for $\gamma > 1$ the shift of the eigenvalues is vanishingly small, and thus they remain close to their initial configuration $\lambda_i(t=0) = a_i$: the spectral density is then given by $p_a$ as expected. Conversely, for $\gamma < 1$, the $a_i = O(1)$ are negligible with respect to the shift of $O(N^{1-\gamma})$, and thus the eigenvalues essentially follow the evolution of a pure Wishart process, whose stationary distribution is known to be the Marčenko–Pastur distribution [98].

For the eigenvectors, one can perform a similar analysis. At short times, for $n \neq i$, the only non-zero term in equation (95) is the last one:

$$ \frac{d\psi_i(n)}{dt} \approx \sigma \sum_{j \neq i} \delta_{jn} \frac{\zeta_{ij}}{a_i - a_j} = \sigma \frac{\zeta_{in}}{a_i - a_n}, \quad (100) $$

and equation (96) simplifies to $\langle \zeta_{ij}(t)\zeta_{kl}(t') \rangle \approx \delta(t-t') [\lambda_i(t) + \lambda_j(t) - a_i - a_j] (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})$. The variance of the components of the eigenvectors at time $T$ is then given by

$$ \begin{aligned} \langle |\psi_i(n; T)|^2 \rangle &= \left( \frac{\sigma}{a_i - a_n} \right)^2 \int_0^T dt \int_0^T dt' \langle \zeta_{in}(t) \zeta_{in}(t') \rangle \\ &\approx \frac{\sigma^2}{(a_i - a_n)^2} \int_0^T dt (\lambda_i(t) + \lambda_n(t) - a_i - a_n). \end{aligned} \quad (101) $$

Now, using $\lambda_i(t) - a_i \approx \sigma^2 Mt$, and the fact that the mean level spacing scales as $a_i - a_n \sim 1/N$, we obtain that after a time $T \sim N^{-\gamma}$

$$ \langle |\psi_i(n; T)|^2 \rangle \sim N^2 \times MT^2 \sim N^3 T^2 \sim N^{3-2\gamma}. \quad (102) $$

Therefore, the eigenvectors are delocalized whenever $\gamma < 3/2$, and are localized otherwise.
7. Conclusions

In this paper, we introduced a new variant of the RP model, which we called the
Wishart-RP ensemble. This model is defined as the sum of a diagonal matrix and
a Wishart matrix, as shown in equations (2) and (3). We provided a comprehensive
analysis of its spectral and localization properties in the limit of large matrix size,
combining several complementary analytical approaches. In particular, we characterized
the phase diagram of the model using perturbation theory (section 3), the cavity method
 subsection 4), the replica formalism (section 5) and a Dyson Brownian motion approach
 subsection 6).

The WRP ensemble displays two distinct transitions (see figure 1): a transition in
the spectral density at $\gamma = 1$, and a localization transition in the eigenvector statistics
at $\gamma = 3/2$. In the intermediate regime ($1 < \gamma < 3/2$), the eigenvectors are neither fully
localized nor extended, but instead occupy a fractal support of size $N^D$, with fractal
dimension $D = 3 - 2\gamma < 1$, which is thus much smaller than $N$.

We also analyzed the spectral correlations through the full counting statistics and
the level compressibility defined in equation (10). The key result of our work is that, in
the intermediate phase ($1 < \gamma < 3/2$), for energy scales between the mean level spacing
and the Thouless energy ($1/N \ll E_T \sim N^{2-2\gamma} \ll 1$), the level compressibility follows the
same scaling function as in other RP-type models (see equation (13)). This agreement,
observed both analytically and numerically, strongly supports the hypothesis that spec-
tral correlations on this scale are (super)universal. While their behavior at smaller and
larger energy scales depends on the specific form of the underlying matrix distributions,
the intermediate-scale correlations appear to be independent of any microscopic detail.

This finding motivates further numerical investigations of the crossover function (10) in realistic many-body disordered quantum systems exhibiting (multi)fractal phases, to test whether the same universal behavior persists in those contexts as well. The most natural framework to begin this investigation is provided, in our view, by the quantum random energy model [22–25, 53].

Finally, several open questions emerge from our study. First, Wishart matrices often
describe covariance data where the number of samples greatly exceeds the number of
variables [75, 76]; this motivates extending the WRP ensemble to the regime $M \ll N$,
e.g. $M \sim N^{\gamma'}$ with $\gamma' < 1$. Second, both the GOE and Wishart ensembles share the
key feature of having a spectral density of the form $\exp[-\mathrm{Tr}V(M)]$, leading to Haar-
distributed eigenvectors. This common structure may underlie the observed universality
and suggests possible extensions to generic matrix models with potentials $V(M)$. Third,
having found that correlations between the entries of $\mathbf{B}$ preserve the conjectured uni-
versality of the level compressibility, it would be interesting to explore what happens
when correlations are also introduced among the entries of $\mathbf{A}$ [72, 73]. Another direction
is to study spectral correlations in models with explicitly multifractal eigenstates, for
example by considering sums of several random matrices with distinct fractal dimension
spectra. Finally, it would be highly desirable to develop a more systematic framework
to study spectral correlations, possibly using tools from free probability theory, which
would constitute a promising avenue for future research.
Acknowledgments

We thank Ivan Khaymovich, Guilhem Semerjian and Pierpaolo Vivo for useful discussions, and Victor Sitzmann for his help with the numerical calculations. We acknowledge financial support from the ANR research grant ManyBodyNet ANR-24-CE30-5851. G S acknowledges support from ANR Grant No. ANR-23-CE30-0020-01 EDIPS.

Appendix A. Scaling function of the level compressibility

In this appendix, we provide the details of the computation that leads to the closed-form scaling function (13), assumed by the level compressibility $\chi(E)$ in the vicinity of the Thouless energy $E_T$. To this end, we analyze the integral that enters the definition of $\chi(x)$ given in equation (61). We introduce the function

$$F_{\chi}(b) = \int_{-\infty}^{\infty} d\tilde{a} [\operatorname{atan}(b + \tilde{a}) - \operatorname{atan}(\tilde{a} - b)]^2, \qquad (A.1)$$

such that $\chi(x) \approx \frac{cp_a(0)}{2\pi x} F_{\chi}(x/(c\pi p_a(0)))$. We first notice that $F_{\chi}(0) = 0$ and compute its derivative $F'_{\chi}(b)$, which reads

$$F'_{\chi}(b) = 2 \int_{-\infty}^{\infty} d\tilde{a} \left( \frac{1}{1 + (\tilde{a} - b)^2} + \frac{1}{1 + (\tilde{a} + b)^2} \right) (\operatorname{atan}(\tilde{a} + b) - \operatorname{atan}(\tilde{a} - b)). \quad (A.2)$$

By using the symmetry of the interval of integration together with $\operatorname{atan}(-x) = -\operatorname{atan}(x)$, the derivative $F'_{\chi}(b)$ can be simply re-written as

$$F'_{\chi}(b) = 4 \int_{-\infty}^{\infty} du \frac{\operatorname{atan}(u)}{1 + (u - 2b)^2}. \qquad (A.3)$$

To proceed, it is convenient to use the integral representation

$$\frac{1}{1+x^2} = \frac{1}{2} \int_{-\infty}^{\infty} dk \, e^{ikx-|k|}. \qquad (A.4)$$

We then insert this integral representation into equation (A.3) and use the identity

$$\int_{-\infty}^{\infty} du \operatorname{atan}(u) e^{iku} = \frac{i\pi}{k} e^{-|k|} \qquad (A.5)$$

(which can be shown using integration by parts), to compute the integral over $u$ in equation (A.3). This leads to

$$F'_{\chi}(b) = 2\pi \int_{-\infty}^{\infty} \frac{dk}{k} \sin(2kb) e^{-2|k|} = 4\pi \operatorname{atan}(b), \qquad (A.6)$$

which can be easily proven by taking a derivative with respect to $b$. Finally, $F_{\chi}(b)$ can be obtained by integrating (A.6) and using $F_{\chi}(0) = 0$, yielding
$$F_{\chi}(b) = \int_{0}^{b} dx \operatorname{atan}(x) = b \operatorname{atan}(b) - \frac{1}{2} \ln(1 + b^2). \qquad (A.7)$$

Using $\chi(x) \approx \frac{c p_a(0)}{2\pi x} F_{\chi}(x/(c\pi p_a(0)))$, this finally leads to equation (13).

# Appendix B. Details of the replica calculation

Here we provide details of the derivations presented in section 5.

## B.1. Density of states

In this appendix, we provide details on the derivation of the DoS using the replica formalism. We start by calculating the average of the partition function replicated $n$ times, according to equations (66) and (67):

$$ \begin{aligned} \langle \mathcal{Z}^n(\lambda) \rangle &= \left\langle \int \prod_{\alpha=1}^{n} d\boldsymbol{r}_{\alpha} e^{-\frac{i}{2} \sum_{\alpha=1}^{n} \boldsymbol{r}_{\alpha}^{T}(\lambda \mathbf{1} - \mathbf{H})\boldsymbol{r}_{\alpha}} \right\rangle_{\mathbf{H}} \\ &= \int \prod_{\alpha=1}^{n} d\boldsymbol{r}_{\alpha} e^{-\frac{i}{2}\lambda \sum_{i=1}^{N} \sum_{\alpha=1}^{n} r_{i\alpha}^2} \left\langle e^{\frac{i}{2} \sum_{i,j=1}^{N} \sum_{\alpha=1}^{n} r_{i\alpha} r_{j\alpha} H_{ij} r_{j\alpha}} \right\rangle_{\mathbf{H}}. \end{aligned} \qquad (B.1) $$

Then we insert the definition of the model $H_{ij} = a_i \delta_{ij} + \nu M^{-\gamma} \sum_{k=1}^{M} W_{ik} W_{jk}$ (see section 2), where $W_{ij} \sim \mathcal{N}(0,1)$, and $a_i$ are i.i.d. random variables sampled from $p_a(a)$. Using that the matrices $\mathbf{A}$ and $\mathbf{W}$ are independent, we now rewrite the term in brackets as

$$ \begin{aligned} \left\langle e^{\frac{i}{2} \sum_{i,j=1}^{N} \sum_{\alpha=1}^{n} r_{i\alpha} H_{ij} r_{j\alpha}} \right\rangle_{\mathbf{H}} &= \left\langle e^{\frac{i}{2} \sum_{i=1}^{N} \sum_{\alpha=1}^{n} r_{i\alpha}^2 a_i} \right\rangle_{\mathbf{A}} \left\langle e^{\frac{i\nu}{2M^{\gamma}} \sum_{i,j=1}^{N} \sum_{\alpha=1}^{n} r_{i\alpha} \sum_{k=1}^{M} W_{ik} W_{jk} r_{j\alpha}} \right\rangle_{\mathbf{W}} \\ &= \left\langle e^{\frac{i}{2} \sum_{i=1}^{N} \sum_{\alpha=1}^{n} r_{i\alpha}^2 a_i} \right\rangle_{\mathbf{A}} \left\langle e^{\frac{i\nu}{2M^{\gamma}} \sum_{k=1}^{M} \sum_{\alpha=1}^{n} \left(\sum_{i=1}^{N} r_{i\alpha} W_{ik}\right)^2} \right\rangle_{\mathbf{W}}. \end{aligned} \qquad (B.2) $$

Given that the measure is Gaussian, the average over $\mathbf{W}$ could in principle be computed immediately (see e.g. [99], where this strategy is applied to the case of Wishart product matrices). However, this gives rise to determinants, which are less convenient in our context in view of the following calculation. An alternative strategy, applied for instance in [100] to the case of diluted Wishart matrices, is to introduce $M \times n$ Hubbard–Stratonovich transformations, with auxiliary variables $u_{k\alpha}$, to decouple the squared term before taking the average. This gives (omitting the $(1/2\pi)^{nM/2}$ prefactor)

$$ \left\langle e^{\frac{i\nu}{2M^{\gamma}} \sum_{k=1}^{M} \sum_{\alpha=1}^{n} \left(\sum_{i=1}^{N} r_{i\alpha} W_{ik}\right)^2} \right\rangle_{\mathbf{W}} \propto \int \prod_{k=1}^{M} \prod_{\alpha=1}^{n} du_{k\alpha} e^{-\frac{u_{k\alpha}^2}{2}} \left\langle e^{\frac{\sqrt{\nu}}{\sqrt{1M^{\gamma}}} \sum_{k=1}^{M} \sum_{\alpha=1}^{n} u_{k\alpha} \sum_{i=1}^{N} r_{i\alpha} W_{ik}} \right\rangle_{\mathbf{W}}, \qquad (B.3) $$
where the average over the Gaussian variables $W_{ik}$ can be easily computed as

$$ \left\langle e^{\frac{\sqrt{\nu}}{\sqrt{iM^{\gamma}}} \sum_{k=1}^{M} \sum_{\alpha=1}^{n} u_{k\alpha} \sum_{i=1}^{N} r_{i\alpha} W_{ik}} \right\rangle_{\mathbf{W}} = e^{\frac{\nu}{2iM^{\gamma}} \sum_{i=1}^{N} \sum_{k=1}^{M} \left( \sum_{\alpha=1}^{n} u_{k\alpha} r_{i\alpha} \right)^2}. \quad (B.4) $$

We now introduce the normalized densities $\phi(\vec{u})$ and $\psi(\vec{r})$, defined as

$$ \phi(\vec{u}) = \frac{1}{M} \sum_{k=1}^{M} \prod_{\alpha=1}^{n} \delta(u_{\alpha} - u_{k\alpha}) = \frac{1}{M} \sum_{k=1}^{M} \delta(\vec{u} - \vec{u}_k), \quad (B.5) $$

$$ \psi(\vec{r}) = \frac{1}{N} \sum_{i=1}^{N} \prod_{\alpha=1}^{n} \delta(r_{\alpha} - r_{i\alpha}) = \frac{1}{N} \sum_{i=1}^{N} \delta(\vec{r} - \vec{r}_i), \quad (B.6) $$

where $\vec{u}_k$ and $\vec{r}_i$ are the $n$-dimensional vectors extracted from the lines of $u_{k\alpha}$ and $r_{i\alpha}$, i.e. $\vec{u}_k = (u_{k\alpha})_{\alpha=1}^n$ and $\vec{r}_i = (r_{i\alpha})_{\alpha=1}^n$. (By contrast, in equation (B.1) we have used boldface to denote a vector $\mathbf{r} \in \mathbb{R}^N$.) Inserting these expressions, and performing simplifications similar to the ones applied in [61] to the case of the GRP model, we finally recover the path-integral representation reported in equations (68) and (69).

We can now evaluate the average of the replicated partition function through the saddle-point method in the limit $N, M \to \infty$ but with the ratio $c = N/M$ kept fixed. The condition of minimization of the action, namely

$$ \frac{\delta S_n}{\delta \psi} = \frac{\delta S_n}{\delta \phi} = \frac{\delta S_n}{\delta \hat{\psi}} = \frac{\delta S_n}{\delta \hat{\phi}} = 0, \quad (B.7) $$

yields the following four saddle-point equations:

$$ \hat{\psi}(\vec{r}) = -\frac{1}{2}\eta \int d\vec{u} \, \phi(\vec{u}) (\vec{u} \cdot \vec{r})^2, \quad (B.8) $$

$$ \psi(\vec{r}) = \frac{\varphi_a(-\vec{r}^2/2) \exp[-\frac{i}{2}\lambda\vec{r}^2 + i\hat{\psi}(\vec{r})]}{\int d\vec{r}' \varphi_a(-\vec{r}'^2/2) \exp[-\frac{i}{2}\lambda\vec{r}'^2 + i\hat{\psi}(\vec{r}')]}, \quad (B.9) $$

$$ \hat{\phi}(\vec{u}) = -\frac{1}{2}c\eta \int d\vec{r} \, \psi(\vec{r}) (\vec{u} \cdot \vec{r})^2, \quad (B.10) $$

$$ \phi(\vec{u}) = \frac{\exp[-\frac{1}{2}\vec{u}^2 + i\hat{\phi}(\vec{u})]}{\int d\vec{u}' \exp[-\frac{1}{2}\vec{u}'^2 + i\hat{\phi}(\vec{u}')]}, \quad (B.11) $$

where we recall the definition of $\eta = \nu M^{1-\gamma}$. Now we note that, according to the Edwards–Jones formula (65), the spectral density can be recovered as

$$ \rho(\lambda) = -\frac{2}{N\pi} \lim_{\epsilon \to 0^+} \mathrm{Im} \frac{\partial}{\partial \lambda} \lim_{n \to 0} \frac{1}{n} \ln \langle \mathcal{Z}^n(\lambda) \rangle. \quad (B.12) $$
Using the saddle-point construction and equation (68), we can express $\langle Z^n(\lambda) \rangle$ as

$$ \langle Z^n(\lambda) \rangle \approx \exp \left\{ -\sqrt{NM} S_n \left[ \phi^*, \hat{\phi}^*, \psi^*, \hat{\psi}^*; \lambda \right] \right\}, \quad (B.13) $$

where the superscript ‘star’ (introduced here for clarity, but omitted hereafter) indicates that these fields are solutions of the saddle-point equations (B.8)–(B.11). In particular, we note that there is only one term in the action (69) in which $\lambda$ appears explicitly, and it only involves the field $\hat{\psi}$. This implies that, in order to derive the spectral density, we only need to find the solution for $\hat{\psi}$—indeed, contributions coming from the implicit derivatives of the other fields with respect to $\lambda$ (calculated according the standard chain rule) vanish by construction at the saddle point due to equation (B.7) (see [61]). Thus, our strategy will be to find a self-consistent equation for $\hat{\psi}$, which can be achieved as follows:

(i) First, we insert equation (B.10) into equation (B.11) to eliminate $\hat{\phi}$ and get an expression for $\phi$ as a functional of $\psi$:

$$ \phi(\vec{u}) = \frac{1}{z_{\phi}} \exp \left[ -\frac{1}{2} \vec{u}^2 - i \frac{c\eta}{2} \int d\vec{r} \, \psi(\vec{r}) (\vec{u} \cdot \vec{r})^2 \right], \quad (B.14) $$

where $z_{\phi} = \int d\vec{u} \exp \left[ -\frac{1}{2}\vec{u}^2 - i \frac{c\eta}{2} \int d\vec{r} \, \psi(\vec{r}) (\vec{u} \cdot \vec{r})^2 \right]$.

(ii) Second, we insert equation (B.9) into equation (B.14) to get an expression for $\phi$ as a functional of $\hat{\psi}$:

$$ \phi(\vec{u}) = \frac{1}{z_{\phi}} \exp \left[ -\frac{1}{2}\vec{u}^2 - i \frac{c\eta}{2z_{\psi}} \int d\vec{r} \, \varphi_a \left(-\frac{\vec{r}^2}{2}\right) e^{-\frac{i}{2}\lambda\vec{r}^2 + i\hat{\psi}(\vec{r})} (\vec{u}' \cdot \vec{r})^2 \right], \quad (B.15) $$

where $z_{\psi} = \int d\vec{r} \, \varphi_a(-\vec{r}^2/2) \exp \left[ -\frac{i}{2}\lambda\vec{r}^2 + i\hat{\psi}(\vec{r}) \right]$.

(iii) Finally, we insert equation (B.15) into equation (B.8) and deduce the self-consistency equation for $\hat{\psi}$:

$$ \hat{\psi}(\vec{r}) = -\frac{\eta}{2z_{\psi}} \int d\vec{u} e^{-\frac{\vec{u}^2}{2} - \frac{i\eta}{2z_{\psi}} \int d\vec{r} \, \varphi_a \left(-\frac{\vec{r}^2}{2}\right) \exp \left[-\frac{i\lambda\vec{r}^2}{2} + i\hat{\psi}(\vec{r})\right] (\vec{u}' \cdot \vec{r})^2} (\vec{u} \cdot \vec{r})^2. \quad (B.16) $$

In the following, we shall try to find a solution of equation (B.16) in the form of a rotationally invariant Ansatz, i.e. $\hat{\psi}(\vec{r}) = \hat{\psi}(r)$, which only depends on the norm of $\vec{r}$ in replica space. Using the identity $\int d\Omega_n (\vec{u} \cdot \vec{r})^2 = \frac{(ur)^2}{n} \int d\Omega_n$, where $d\Omega_n$ is the differential of the n-dimensional solid angle in spherical coordinates, we obtain

$$ \hat{\psi}(r) = -r^2 \frac{\eta}{2nz_{\psi}} \int du \, u^{n+1} e^{-\frac{u^2}{2} - \frac{i\eta}{2nz_{\psi}} \int d\Omega_n u^2 \int dr' r'^{n+1} \varphi_a \left(-\frac{r'^2}{2}\right) \exp \left[-\frac{i\lambda r'^2}{2} + i\hat{\psi}(r')\right]} \quad (B.17) $$

This expression can then be simplified by defining

$$ J(r; \lambda) = \varphi_a \left(-\frac{r^2}{2}\right) \exp \left[ -\frac{i}{2}\lambda r^2 + i\hat{\psi}(r) \right], \quad (B.18) $$
with $J'(r; \lambda) = \partial_r J(r; \lambda)$. Rewriting $z_\psi = \int d\Omega_n \int dr r^{n-1} J(r; \lambda)$ and integrating by parts
$$z_\psi = \frac{\int d\Omega_n}{n} \int dr r^n J'(r; \lambda), \text{ we obtain}$$

$$\hat{\psi}(r) = -r^2 \frac{\eta \int d\Omega_n}{2nz_\phi} \int du \, u^{n+1} \exp \left[ -\frac{1}{2}u^2 + i \frac{c\eta}{2} u^2 \frac{\int dr' r'^{n+1} J(r'; \lambda)}{\int dr' r'^n J'(r'; \lambda)} \right]. \quad (B.19)$$

This expression can be further simplified by defining

$$F_n(u; \lambda) = \exp \left[ -\frac{1}{2}u^2 + i \frac{c\eta}{2} u^2 \frac{\int dr \, r^{n+1} J(r; \lambda)}{\int dr \, r^n J'(r; \lambda)} \right], \quad (B.20)$$

and $F'_n(u; \lambda) = \partial_u F_n(u; \lambda)$. Again, rewriting $z_\phi = \int d\Omega_n \int du \, u^{n-1} F_n(u; \lambda)$ and integrating by parts $z_\phi = \frac{\int d\Omega_n}{n} \int du \, u^n F'_n(u; \lambda)$, we have

$$\hat{\psi}(r) = \frac{\eta}{2} r^2 \frac{\int du \, u^{n+1} F_n(u; \lambda)}{\int du \, u^n F'_n(u; \lambda)}. \quad (B.21)$$

We can now take the limit $n \to 0$ and arrive at

$$\hat{\psi}(r) = \frac{\eta}{2} r^2 \frac{\int du \, u F_0(u; \lambda)}{\int du \, F'_0(u; \lambda)}. \quad (B.22)$$

Let us now go back to the Edwards-Jones formula and the replica trick in equations (65) and (67), which allow us to compute the spectral density as

$$\begin{align*}
\rho(\lambda) &= -\frac{2}{\pi N} \lim_{\epsilon \to 0^+} \mathrm{Im} \frac{\partial}{\partial \lambda} \lim_{n \to 0} \frac{\ln \langle \mathcal{Z}^n(\lambda_\epsilon) \rangle}{n} = \frac{2}{\pi N} \lim_{\epsilon \to 0^+} \mathrm{Im} \lim_{n \to 0} \frac{1}{n} \sqrt{NM} \frac{\partial}{\partial \lambda} S_n^{\mathrm{sp}} \\
&= \frac{2}{\pi \sqrt{c}} \lim_{\epsilon \to 0^+} \mathrm{Im} \lim_{n \to 0} \frac{-i\sqrt{c}}{2} \frac{1}{n} \frac{\int dr \, r^{n+1} J(r; \lambda_\epsilon)}{\int dr \, r^{n-1} J(r; \lambda_\epsilon)} \\
&= \frac{1}{\pi} \lim_{\epsilon \to 0^+} \mathrm{Im} i \frac{\int dr \, r J(r; \lambda_\epsilon)}{\int dr \, J'(r; \lambda_\epsilon)} \equiv \frac{1}{\pi} \lim_{\epsilon \to 0^+} \mathrm{Im} G(\lambda_\epsilon), \quad (B.23)
\end{align*}$$

where we have identified the resolvent

$$G(\lambda) = i \frac{\int dr \, r J(r; \lambda)}{\int dr \, J'(r; \lambda)}. \quad (B.24)$$

To obtain a closed equation for the resolvent, we first rewrite equation (B.20) as

$$F_0(u; \lambda) = \exp \left[ -\frac{1}{2} u^2 + \frac{c\eta}{2} u^2 G(\lambda) \right]. \quad (B.25)$$

From this expression, if we assume that $G(\lambda)$ has a negative real part on the support of $\rho$ (to be checked a posteriori), then we have that $F_0(u \to \infty; \lambda) = 0$ and $F_0(u \to 0; \lambda) = 1$. This allows us to simplify the denominator of equation (B.22), and then by applying
the change of variables $y = u^2/2$, we find

$$ 
\begin{aligned}
\hat{\psi}(r) &= \frac{\eta}{2} r^2 \int_0^\infty du \, u \exp \left[ -\frac{1}{2} u^2 + \frac{c\eta}{2} u^2 G(\lambda) \right] \\
&= \frac{\eta}{2} r^2 \int_0^\infty dy \, \exp \left[ -y \left( 1 - c\eta G(\lambda) \right) \right] = \frac{\eta}{2} r^2 \frac{1}{1 - c\eta G(\lambda)}.
\end{aligned} 
\quad (B.26) $$

Next, since $G(\lambda)$ has a positive imaginary part, from equations (B.18) and (B.26) we deduce that $J(r \to \infty; \lambda) = 0$ and $J(r \to 0; \lambda) = 1$. These limits let us simplify the denominator of equation (B.24), so that by applying the change of variables $z = r^2/2$ and inserting equation (B.18), we get

$$ 
\begin{aligned}
G(\lambda) &= i \int_0^\infty dr \, r \, \varphi_a \left(-\frac{r^2}{2}\right) \exp \left[ -\frac{i}{2} \lambda r^2 + i \hat{\psi}(r) \right] \\
&= i \int_0^\infty dr \, r \, \varphi_a \left(-\frac{r^2}{2}\right) \exp \left[ -\frac{i}{2} \lambda r^2 + \frac{i}{2} r^2 \frac{\eta}{1 - c \eta G(\lambda)} \right].
\end{aligned} 
\quad (B.27) $$

Finally, upon changing variables as $z = r^2/2$, we recover equation (72).

## B.2. Full counting statistics and level compressibility

To obtain the action (82), we follow very similar steps to the ones presented in the previous section for the DoS: we start from equation (79), we perform the Hubbard-Stratonovich transformation, we average over the Gaussian measure and we insert the fields as in equations (B.8)–(B.11) to obtain equation (81).

Now, to evaluate the action (82) at the saddle point, we first write the saddle-point equations for the fields, which read

$$ \hat{\psi}(\vec{r}) = -\frac{1}{2}\eta \int d\vec{u} \, \phi(\vec{u}) M(\vec{u}, \vec{r}), \quad (B.28) $$

$$ \hat{\phi}(\vec{u}) = -\frac{1}{2}c\eta \int d\vec{r} \, \psi(\vec{r}) M(\vec{u}, \vec{r}), \quad (B.29) $$

$$ \psi(\vec{r}) = \frac{1}{Z_{\psi}} \varphi_a \left( -\frac{1}{2} \vec{r} \hat{L} \vec{r} \right) \exp \left[ -\frac{i}{2} \vec{r} \hat{\Lambda} \vec{r} + i \hat{\psi}(\vec{r}) \right], \quad (B.30) $$

$$ \phi(\vec{u}) = \frac{1}{Z_{\phi}} \exp \left[ -\frac{1}{2} \vec{u}^2 + i \hat{\phi}(\vec{u}) \right], \quad (B.31) $$

where we introduced $M(\vec{u}, \vec{r}) \equiv (\vec{u}\hat{L}\vec{r})^2$. Inserting the first two equations back into equation (82), the action at the saddle point simplifies to

$$ S_{n\pm} [\phi, \psi; \hat{\Lambda}] = -\sqrt{c} \ln Z_{\psi} - \frac{1}{\sqrt{c}} \ln Z_{\phi} - \frac{i}{2} \sqrt{c\eta} \int d\vec{r} d\vec{u} \phi(\vec{u}) \psi(\vec{r}) M(\vec{u}, \vec{r}), \quad (B.32) $$
where

$$ 
\begin{aligned}
Z_{\psi} &= \int d\vec{r} \, \varphi_a \left(-\frac{1}{2}\vec{r}\hat{L}\vec{r}\right) \exp \left[-\frac{i}{2}\vec{r}\hat{\Lambda}\vec{r} + i\hat{\psi}(\vec{r})\right], \\
Z_{\phi} &= \int d\vec{u} \, \exp \left[-\frac{1}{2}\vec{u}^2 + i\hat{\phi}(\vec{u})\right].
\end{aligned} 
\quad (B.33) $$

On the other hand, by inserting the first two saddle-point equations into the last two, one can eliminate the dependence on the conjugated fields $\hat{\psi}$ and $\hat{\phi}$, and obtain a set of two equations for $\psi$ and $\phi$ only:

$$ \psi(\vec{r}) = \frac{1}{Z_{\psi}} \varphi_a \left(-\frac{1}{2}\vec{r}\hat{L}\vec{r}\right) \exp \left[-\frac{i}{2}\vec{r}\hat{\Lambda}\vec{r} - \frac{i}{2}\eta \int d\vec{u} \, \phi(\vec{u}) M(\vec{u}, \vec{r})\right], \quad (B.34) $$

$$ \phi(\vec{u}) = \frac{1}{Z_{\phi}} \exp \left[ -\frac{1}{2}\vec{u}^2 - \frac{i}{2}c\eta \int d\vec{r} \, \psi(\vec{r}) M(\vec{u}, \vec{r}) \right]. \quad (B.35) $$

To make progress, we now introduce an $n \times n$ matrix $\hat{K}$, such that

$$ \vec{u}\hat{K}\vec{u} = \int d\vec{r} \, \psi(\vec{r}) M(\vec{u}, \vec{r}), \quad (B.36) $$

which allows us to rewrite equations (B.34) and (B.35) as

$$ \phi(\vec{u}) = \frac{1}{Z_{\phi}} \exp \left[ -\frac{1}{2}\vec{u} \left( \hat{1} + iq\eta\hat{K} \right) \vec{u} \right], \quad (B.37) $$

$$ \psi(\vec{r}) = \frac{1}{Z_{\psi}} \varphi_a \left( -\frac{1}{2}\vec{r}\hat{L}\vec{r} \right) \exp \left\{ -\frac{i}{2}\vec{r} \left[ \hat{\Lambda} + \eta\hat{L} \left( \hat{1} + iq\eta\hat{K} \right)^{-1} \hat{L} \right] \vec{r} \right\}. \quad (B.38) $$

Inserting equation (B.38) into equation (B.36) then gives the following self-consistent equation for $\hat{K}$:

$$ \vec{u}\hat{K}\vec{u} = \frac{1}{Z_{\psi}} \int d\vec{r} \, \varphi_a \left( -\frac{1}{2}\vec{r}\hat{L}\vec{r} \right) e^{-\frac{i}{2}\vec{r}[\hat{\Lambda}+\eta\hat{L}(\hat{1}+iq\eta\hat{K})^{-1}\hat{L}]\vec{r}} M(\vec{u}, \vec{r}), \quad (B.39) $$

which is actually equivalent to equation (85) after introducing the auxiliary matrix $\hat{C}$. Using equations (B.37) and (B.38), we can now express $Z_{\phi}$ and $Z_{\psi}$ in equation (B.33) as

$$ Z_{\phi} = (2\pi)^{n/2} \times \left[ \det (\hat{1} + iq\eta\hat{K}) \right]^{-1/2}, \quad (B.40) $$

$$ Z_{\psi} = \int_{-\infty}^{\infty} da \, p_a(a) \left[ \det (\hat{C}^{-1} - ia\hat{L}) \right]^{-1/2}. \quad (B.41) $$

Furthermore, the interaction term in the action (B.32) can be rewritten as
The Wishart–Rosenzweig–Porter random matrix ensemble

$$ \begin{aligned} \int d\vec{r} d\vec{u} \phi(\vec{u}) \psi(\vec{r}) M(\vec{u}, \vec{r}) &= \int d\vec{u} \phi(\vec{u}) \int d\vec{r} \psi(\vec{r}) M(\vec{u}, \vec{r}) = \int d\vec{u} \phi(\vec{u}) \vec{u} \hat{K} \vec{u} \\ &= \frac{1}{Z_{\phi}} \int d\vec{u} \exp \left[ -\frac{1}{2} \vec{u} (\hat{1} + i q \eta \hat{K}) \vec{u} \right] \vec{u} \hat{K} \vec{u} \\ &= \mathrm{Tr} \left[ \hat{K} (\hat{1} + i q \eta \hat{K})^{-1} \right]. \end{aligned} \quad (B.42) $$

Combining these results, the expression in equation (B.32) finally simplifies to the one in equation (84).

# Appendix C. Details of the Dyson Brownian motion calculation

Here we provide details of the derivations presented in section 6.

## C.1. First-order perturbation theory for the eigenvalues

We begin by detailing the steps leading from equation (93) to the evolution equation (94) for the eigenvalues. To this end, we resort to perturbation theory. At first order, the latter tells us that the perturbation $\delta\mathbf{H}(t)$ shifts the eigenvalues of $\mathbf{H}(t)$ by $\delta^{(1)}\lambda_i = \langle\psi_i(t)|\delta\mathbf{H}(t)|\psi_i(t)\rangle$, where $|\psi_i(t)\rangle$ is the $i$th eigenvector of $\mathbf{H}(t)$ at time $t$.

First, let us study the effect of the term $\mathbf{g}(t)\mathbf{X}(t)^T$ in $\delta\mathbf{H}(t)$ on the eigenvalues:

$$ \langle \psi_i(t) | \mathbf{g}(t) \mathbf{X}(t)^T | \psi_i(t) \rangle = \sum_{n,m} \sum_l \psi_i(n;t) g_{nl}(t) X_{ml}(t) \psi_i(m;t) \equiv \xi_i(t), \quad (C.1) $$

where we introduced the random variable $\xi_i(t)$. As this random variable is itself a sum of random variables of finite variance, by virtue of the central limit theorem we expect it to be Gaussian distributed in the large-$N$ limit. Adopting the Itô convention, its average is given by

$$ \begin{aligned} \langle \xi_i(t) \rangle &= \left\langle \sum_{n,m} \sum_l \psi_i(n;t) g_{nl}(t) X_{ml}(t) \psi_i(m;t) \right\rangle \\ &= \sum_{n,m} \sum_l \psi_i(n;t) \langle g_{nl}(t) \rangle X_{ml}(t) \psi_i(m;t) = 0, \end{aligned} \quad (C.2) $$

while its variance reads

$$ \begin{aligned} \langle \xi_i(t) \xi_j(t') \rangle &= \left\langle \sum_{n,m,n',m'} \sum_{l,l'} \psi_i(n;t) g_{nl}(t) X_{ml}(t) \psi_i(m;t) \psi_i(n';t') g_{n'l'}(t') X_{m'l'}(t') \psi_i(m';t') \right\rangle \\ &= \sum_{n,m,n',m'} \sum_{l,l'} \psi_i(n;t) \psi_i(m;t) \psi_i(n';t') \psi_i(m';t') \langle g_{nl}(t) g_{n'l'}(t') \rangle X_{ml}(t) X_{m'l'}(t') \\ &= \sigma^2 \mathrm{d}t \, \delta(t-t') \sum_{n,m,n',m'} \sum_{l,l'} \psi_i(n;t) \psi_i(m;t) \psi_i(n';t') \psi_i(m';t') X_{ml}(t) X_{m'l'}(t') \delta_{nl} \delta_{n'l'} \\ &= \sigma^2 \mathrm{d}t \, \delta(t-t') \sum_{n,m,m'} \sum_l \psi_i(n;t) \psi_i(m)(t) \psi_i(n;t) \psi_i(m';t) X_{ml}(t) X_{m'l}(t). \end{aligned} \quad (C.3) $$
We now note that, because of the normalization of the eigenvectors $\sum_n \psi_j(n;t)\psi_i(n;t) = \langle i(t)|j(t)\rangle = \delta_{ij}$, the variance can be further simplified as

$$ 
\begin{aligned}
\langle \xi_i(t) \xi_j(t') \rangle &= \sigma^2 dt \delta_{ij} \delta(t-t') \sum_{m,m'} \sum_l \psi_i(m;t) \psi_i(m';t) X_{ml}(t) X_{m'l}(t) \\
&= \sigma^2 dt \delta_{ij} \delta(t-t') \langle \psi_i(t) | \mathbf{X}(t) \mathbf{X}(t)^T | \psi_i(t) \rangle \\
&= \sigma^2 dt \delta_{ij} \delta(t-t') \langle \psi_i(t) | \mathbf{H}(t) - \mathbf{A} | \psi_i(t) \rangle \\
&= \sigma^2 dt \delta_{ij} \delta(t-t') (\lambda_i(t) - \langle a \rangle_i(t)),
\end{aligned}
\quad (C.4) $$

where $\langle a \rangle_i(t)$ was given in equation (97).

Next, following similar steps as in equation (C.1), one can show that the contribution of $\mathbf{X}(t)\mathbf{g}(t)^T$ is the same as that of $\mathbf{g}(t)\mathbf{X}(t)^T$. Finally, the contribution of the term $\mathbf{g}(t)\mathbf{g}(t)^T$ yields

$$ 
\begin{aligned}
\langle \psi_i(t) | \mathbf{g}(t) \mathbf{g}(t)^T | \psi_i(t) \rangle &= \sum_{n,m} \sum_l \psi_i(n;t) g_{nl}(t) g_{ml}(t) \psi_i(m;t) \\
&= M \sigma^2 dt + \mathcal{O}\left(\left(d t^{3/2}\right)\right).
\end{aligned}
\quad (C.5) $$

## C.2. Second-order perturbation theory for the eigenvalues

The second-order correction to the eigenvalues is given by

$$ 
\begin{aligned}
\delta^{(2)} \lambda_i &= \sum_{j \neq i} \frac{|\langle \psi_j(t) | \delta \mathbf{H}(t) | \psi_i(t) \rangle|^2}{\lambda_i - \lambda_j} \\
&= \sum_{j \neq i} \frac{|\langle \psi_j(t) | \mathbf{g}(t) \mathbf{X}(t)^T + \mathbf{X}(t) \mathbf{g}(t)^T | \psi_i(t) \rangle|^2}{\lambda_i - \lambda_j} + \mathcal{O}(\mathrm{d}t^2) \\
&= \sum_{j \neq i} \left[ |\langle \psi_j(t) | \mathbf{g}(t) \mathbf{X}(t)^T | \psi_i(t) \rangle|^2 + |\langle \psi_j(t) | \mathbf{X}(t) \mathbf{g}(t)^T | \psi_i(t) \rangle|^2 \right. \\
& \qquad \left. + 2 \langle \psi_j(t) | \mathbf{g}(t) \mathbf{X}(t)^T | \psi_i(t) \rangle \langle \psi_j(t) | \mathbf{X}(t) \mathbf{g}(t)^T | \psi_i(t) \rangle \right] \frac{1}{\lambda_i - \lambda_j} + \mathcal{O}(\mathrm{d}t^2). 
\end{aligned}
\quad (C.6) $$

In the second line, the second term is the same as the first one upon switching $i \leftrightarrow j$ in the numerator, while the term in the last line will give $\delta_{ij}$ upon averaging, thus giving no contribution to the sum over $j \neq i$. Retaining only the terms up to order $\mathrm{d}t$ and henceforth, unless otherwise stated, omitting the time dependence at $t$, we obtain

$$ 
\begin{aligned}
\delta^{(2)} \lambda_i &= \sum_{j \neq i} \left[ \sum_{n,m,n',m'} \psi_j(n) g_{nl} X_{ml} \psi_i(m) \psi_j(n') g_{n'l'} X_{m'l'} \psi_i(m') + (i \leftrightarrow j) \right] \frac{1}{\lambda_i - \lambda_j} \\
&= \sigma^2 \mathrm{d}t \sum_{j \neq i} \left[ \sum_{n,m,m'} \sum_l \psi_j(n) X_{ml} \psi_i(m) \psi_j(n) X_{m'l} \psi_i(m') + (i \leftrightarrow j) \right] \frac{1}{\lambda_i - \lambda_j}
\end{aligned}
$$
The Wishart–Rosenzweig–Porter random matrix ensemble

$$ 
\begin{aligned}
&= \sigma^2 dt \sum_{j \neq i} \left[ \sum_{m,m'} \sum_l X_{ml} \psi_i(m) X_{m'l} \psi_i(m') + (i \leftrightarrow j) \right] \frac{1}{\lambda_i - \lambda_j} \\
&= \sigma^2 dt \sum_{j \neq i} \frac{\langle \psi_i | \mathbf{X}\mathbf{X}^T | \psi_i \rangle + \langle \psi_j | \mathbf{X}\mathbf{X}^T | \psi_j \rangle}{\lambda_i - \lambda_j} = \sigma^2 dt \sum_{j \neq i} \frac{\langle \psi_i | \mathbf{H} - \mathbf{A} | \psi_i \rangle + \langle \psi_j | \mathbf{H} - \mathbf{A} | \psi_j \rangle}{\lambda_i - \lambda_j} \\
&= \sigma^2 dt \sum_{j \neq i} \frac{\lambda_i + \lambda_j - \langle a \rangle_i - \langle a \rangle_j}{\lambda_i - \lambda_j}. \qquad (C.7)
\end{aligned}
$$

Gathering all the contributions in equations (C.1), (C.5) and (C.7), we finally get equation (94).

### C.3. First-order perturbation theory for the eigenvectors

We now detail the steps leading from equation (93) to the evolution equation (95) for the eigenvectors. The correction to the nth component of the ith eigenvector at first order reads

$$ \delta^{(1)} \psi_i(n) = \sum_{j \neq i} \psi_j(n) \frac{\langle \psi_j | \mathbf{g}\mathbf{X}^T + \mathbf{X}\mathbf{g}^T + \mathbf{g}\mathbf{g}^T | \psi_i \rangle}{\lambda_i - \lambda_j}. \qquad (C.8) $$

Let us then focus on the first two terms,

$$ \sum_{j \neq i} \psi_j(n) \frac{\langle \psi_j | \mathbf{g}\mathbf{X}^T + \mathbf{X}\mathbf{g}^T | \psi_i \rangle}{\lambda_i - \lambda_j} = \sum_{j \neq i} \psi_j(n) \frac{\xi_{ij}}{\lambda_i - \lambda_j}, \qquad (C.9) $$

where we defined the random variable $\xi_{ij} \equiv \langle \psi_j | \mathbf{g}\mathbf{X}^T + \mathbf{X}\mathbf{g}^T | \psi_i \rangle$. Again, since this is a sum of random variables with finite variance, we expect that in the large-$N$ limit it will converge to a Gaussian random variable of mean value

$$ \langle \xi_{ij}(t) \rangle = \sum_{n,m} \sum_{l} \psi_j(n; t) \langle g_{nl}(t) \rangle X_{ml}(t) \psi_i(m; t) + (i \leftrightarrow j) = 0, \qquad (C.10) $$

and variance

$$ 
\begin{aligned}
\langle \xi_{ij}(t) \xi_{kl}(t') \rangle &= \langle \langle \psi_j | \mathbf{g}\mathbf{X}^T | \psi_i \rangle \langle \psi_l | \mathbf{g}\mathbf{X}^T | \psi_k \rangle + \langle \psi_j | \mathbf{g}\mathbf{X}^T | \psi_i \rangle \langle \psi_l | \mathbf{X}\mathbf{g}^T | \psi_k \rangle + (i \leftrightarrow j) \rangle \\
&= \sigma^2 dt \delta(t-t') [\delta_{jl}(\lambda_i \delta_{ik} - \langle a \rangle_{ik}) + \delta_{jk}(\lambda_i \delta_{il} - \langle a \rangle_{il}) + (i \leftrightarrow j)] \\
&= \sigma^2 dt \delta(t-t') [(\lambda_i + \lambda_j)(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) - (\delta_{ik}\langle a \rangle_{jl} + \delta_{il}\langle a \rangle_{jk} + \delta_{jk}\langle a \rangle_{il} + \delta_{jl}\langle a \rangle_{ik})],
\end{aligned}
\qquad (C.11) $$

where $\langle a \rangle_{ik}(t)$ was given in equation (98).
We finally consider the contribution of $\langle \psi_j | gg^T | \psi_i \rangle$, finding

$$ 
\begin{aligned}
\sum_{j \neq i} \psi_j(n) \frac{\langle \psi_j | gg^T | \psi_i \rangle}{\lambda_i - \lambda_j} &= \sum_{j \neq i} \psi_j(n) \frac{\sum_{m,m'} \sum_l \psi_j(m) \psi_i(m') g_{ml} g_{m'l}}{\lambda_i - \lambda_j} \\
&= M \sigma^2 dt \sum_{j \neq i} \psi_j(n) \frac{\sum_m \psi_j(m) \psi_i(m)}{\lambda_i - \lambda_j} + \mathcal{O}(dt^{3/2}) \\
&= M \sigma^2 dt \sum_{j \neq i} \psi_j(n) \frac{\langle \psi_i | \psi_j \rangle}{\lambda_i - \lambda_j} = 0.. \qquad (C.12)
\end{aligned}
$$

### C.4. Second-order perturbation theory for the eigenvectors

The second-order calculation for the eigenvectors has three contributions:

$$ 
\begin{aligned}
\delta^{(2)} \psi_i(n) &= \sum_{k \neq i} \sum_{l \neq i} \psi_k(n) \frac{\langle \psi_k | \delta \mathbf{H}(t) | \psi_l \rangle \langle \psi_l | \delta \mathbf{H}(t) | \psi_i \rangle}{(\lambda_i - \lambda_k)(\lambda_i - \lambda_l)} \\
&\quad - \sum_{k \neq i} \psi_k(n) \frac{\langle \psi_k | \delta \mathbf{H}(t) | \psi_i \rangle \langle \psi_i | \delta \mathbf{H}(t) | \psi_i \rangle}{(\lambda_i - \lambda_k)^2} \\
&\quad - \frac{1}{2} \psi_i(n) \sum_{k \neq i} \frac{|\langle \psi_k | \delta \mathbf{H}(t) | \psi_i \rangle|^2}{(\lambda_i - \lambda_k)^2}.
\end{aligned}
\qquad (C.13)
$$

The first one is given by

$$ 
\begin{aligned}
& \sum_{k \neq i} \sum_{l \neq i} \psi_k(n) \frac{\langle \psi_k | \mathbf{gX}^T + \mathbf{Xg}^T | \psi_l \rangle \langle \psi_l | \mathbf{gX}^T + \mathbf{Xg}^T | \psi_i \rangle}{(\lambda_i - \lambda_k)(\lambda_i - \lambda_l)} + \mathcal{O}(dt^{3/2}) \\
&= \sigma^2 dt \sum_{\substack{k \neq i \\ l \neq i}} \frac{\delta_{kl} \langle \psi_l | \mathbf{X} \mathbf{X}^T | \psi_i \rangle + \delta_{ik} \langle \psi_l | \mathbf{X} \mathbf{X}^T | \psi_l \rangle + \delta_{il} \langle \psi_k | \mathbf{X} \mathbf{X}^T | \psi_l \rangle + \delta_{ll} \langle \psi_k | \mathbf{X} \mathbf{X}^T | \psi_i \rangle}{[\psi_k(n)]^{-1} (\lambda_i - \lambda_k)(\lambda_i - \lambda_l)} \\
&= \sigma^2 dt \sum_{k \neq i} \psi_k(n) \frac{\lambda_i \delta_{ik} - \langle a \rangle_{ik}}{(\lambda_i - \lambda_k)^2} + \sigma^2 dt \sum_{k \neq i} \sum_{l \neq i} \psi_k(n) \frac{\lambda_i \delta_{ik} - \langle a \rangle_{ik}}{(\lambda_i - \lambda_k)(\lambda_i - \lambda_l)} \\
&= -\sigma^2 dt \sum_{k \neq i} \psi_k(n) \frac{\langle a \rangle_{ik}}{(\lambda_i - \lambda_k)^2} - \sigma^2 dt \sum_{k \neq i} \psi_k(n) \frac{\langle a \rangle_{ik}}{(\lambda_i - \lambda_k)} \sum_{l \neq i} \frac{1}{(\lambda_i - \lambda_l)}.
\end{aligned}
\qquad (C.14)
$$

The second contribution reads

$$ 
\begin{aligned}
& - \sum_{k \neq i} \psi_k(n) \frac{\langle \psi_k | \mathbf{gX}^T + \mathbf{Xg}^T | \psi_i \rangle \langle \psi_i | \mathbf{gX}^T + \mathbf{Xg}^T | \psi_i \rangle}{(\lambda_i - \lambda_k)^2} + \mathcal{O}(dt^{3/2}) \\
&= -\sigma^2 dt \sum_{k \neq i} \psi_k(n) \frac{2\delta_{ik} \langle \psi_i | \mathbf{X} \mathbf{X}^T | \psi_i \rangle + 2\delta_{il} \langle \psi_k | \mathbf{X} \mathbf{X}^T | \psi_i \rangle}{(\lambda_i - \lambda_k)^2} \\
&= -2\sigma^2 dt \sum_{k \neq i} \psi_k(n) \frac{\lambda_i \delta_{ik} - \langle a \rangle_{ik}}{(\lambda_i - \lambda_k)^2} = 2\sigma^2 dt \sum_{k \neq i} \psi_k(n) \frac{\langle a \rangle_{ik}}{(\lambda_i - \lambda_k)^2}.
\end{aligned}
\qquad (C.15)
$$
Finally, the last term in equation (C.13) can be simplified using

$$ \sum_{k \neq i} \frac{|\langle \psi_k | \mathbf{g} \mathbf{X}^T + \mathbf{X} \mathbf{g}^T | \psi_i \rangle|^2}{(\lambda_i - \lambda_k)^2} = \sum_{k \neq i} \frac{\lambda_i + \lambda_k - \langle a \rangle_i - \langle a \rangle_k}{(\lambda_i - \lambda_k)^2}. \qquad (C.16) $$

Gathering all the contributions in equations (C.9) and (C.14)–(C.16), we finally obtain equation (95). The noise variance in equation (96) follows instead from equation (C.11).

## ORCID iDs

Victor Delapalme ID 0009-0002-0492-7965

Leticia F Cugliandolo ID 0000-0002-4986-8164

Grégory Schehr ID 0000-0001-6648-5213

Marco Tarzia ID 0000-0001-8113-7049

Davide Venturelli ID 0000-0001-6952-9568

## References

[1] Livan G, Novaes M and Vivo P 2018 Introduction to Random Matrices (Springer International Publishing) (https://doi.org/10.1007/2F978-3-319-70885-0)

[2] Potters M and Bouchaud J P 2020 A First Course in Random Matrix Theory: For Physicists, Engineers and Data Scientists (Cambridge University Press)

[3] Mingo J A and Speicher R 2017 Free Probability and Random Matrices (Springer)

[4] Srednicki M 1994 Phys. Rev. E 50 888-901

[5] Rigol M, Dunjko V and Olshanii M 2008 Nature 452 854-8

[6] Gorny I V, Mirlin A D and Polyakov D G 2005 Phys. Rev. Lett. 95 206603

[7] Basko D M, Aleiner I L and Altshuler B L 2006 Ann. Phys., NY 321 1126-205

[8] Altman E and Vosk R 2015 Annu. Rev. Condens. Matter Phys. 6 383-409

[9] Nandkishore R and Huse D A 2015 Annu. Rev. Condens. Matter Phys. 6 15-38

[10] Abanin D A and Papić Z 2017 Ann. Phys., Lpz. 529 1700169

[11] Alet F and Laflorencie N 2018 C. R. Physique 19 498-525

[12] Abanin D A, Altman E, Bloch I and Serbyn M 2019 Rev. Mod. Phys. 91 021001

[13] Sierant P, Lewenstein M, Scardicchio A, Vidmar L and Zakrzewski J 2025 Rep. Prog. Phys. 88 026502

[14] Macé N, Alet F and Laflorencie N 2019 Phys. Rev. Lett. 123 180601

[15] De Tomasi G, Khaymovich I M, Pollmann F and Warzel S 2021 Phys. Rev. B 104 024202

[16] Gorny I, Mirlin A, Polyakov D and Burin A 2017 Ann. Phys., Lpz. 529 1600360

[17] Tarzia M 2020 Phys. Rev. B 102 014208

[18] Luitz D J, Laflorencie N and Alet F 2015 Phys. Rev. B 91 081103

[19] Serbyn M, Papić Z and Abanin D A 2017 Phys. Rev. B 96 104201

[20] Tikhonov K S and Mirlin A D 2018 Phys. Rev. B 97 214205

[21] Luitz D J, Khaymovich I M and Lev Y B 2020 SciPost Phys. Core 2 006

[22] Faoro L, Feigel’man M V and Ioffe L 2019 Ann. Phys., NY 409 167916

[23] Baldwin C L and Laumann C R 2018 Phys. Rev. B 97 224201

[24] Biroli G, Facoetti D, Schiró M, Tarzia M and Vivo P 2021 Phys. Rev. B 103 014204

[25] Parolini T and Mossi G 2020 arXiv:2007.00315

[26] Kechedzhi K, Smelyanskiy V, McClean J R, Denchev V S, Mohseni M, Isakov S, Boixo S, Altshuler B and Neven H 2018 Leibniz Int. Proc. Inf. 111 9:1-16

[27] Pino M, Kravtsov V E, Altshuler B L and Ioffe L B 2017 Phys. Rev. B 96 214205

[28] Pino M, Ioffe L B and Altshuler B L 2016 Proc. Natl Acad. Sci. USA 113 536-41

[29] Cugliandolo L F, Schehr G, Tarzia M and Venturelli D 2024 Phys. Rev. B 110 174202

[30] Biroli G, Hartmann A K and Tarzia M 2024 Phys. Rev. B 110 014205

[31] Kravtsov V E, Khaymovich I M, Cuevas E and Amini M 2015 New J. Phys. 17 122002
The Wishart–Rosenzweig–Porter random matrix ensemble

[32] von Soosten P and Warzel S 2019 *Lett. Math. Phys.* **109** 905–22

[33] Facoetti D, Vivo P and Biroli G 2016 *Europhys. Lett.* **115** 47003

[34] Truong K and Ossipov A 2016 *Europhys. Lett.* **116** 37002

[35] Bogomolny E and Sieber M 2018 *Phys. Rev. E* **98** 032139

[36] Tomasi G D, Amini M, Bera S, Khaymovich I M and Kravtsov V E 2019 *SciPost Phys.* **6** 14

[37] Amini M 2017 *Europhys. Lett.* **117** 30003

[38] Pino M, Tabanera J and Serna P 2019 *J. Phys. A: Math. Theor.* **52** 475101

[39] Berkovits R 2020 *Phys. Rev. B* **102** 165140

[40] Kravtsov V E, Khaymovich I M, Altshuler B L and Ioffe L B 2020 arXiv:2002.02979

[41] Khaymovich I M, Kravtsov V E, Altshuler B L and Ioffe L B 2020 *Phys. Rev. Res.* **2** 043346

[42] Monthus C 2017 *J. Phys. A: Math. Theor.* **50** 295101

[43] Biroli G and Tarzia M 2021 *Phys. Rev. B* **103** 104205

[44] Buijsman W and Lev Y B 2022 *SciPost Phys.* **12** 082

[45] Khaymovich I M and Kravtsov V E 2021 *SciPost Phys.* **11** 045

[46] Roy S, Khaymovich I M, Das A and Moessner R 2018 *SciPost Phys.* **4** 025

[47] Wang J, Liu X J, Xianlong G and Hu H 2016 *Phys. Rev. B* **93** 104504

[48] Nosov P A, Khaymovich I M and Kravtsov V E 2019 *Phys. Rev. B* **99** 104203

[49] Duthie A, Roy S and Logan D E 2022 *Phys. Rev. B* **106** L020201

[50] Kutlin A G and Khaymovich I M 2021 *SciPost Phys.* **11** 101

[51] Motamarri V R, Gorsky A S and Khaymovich I M 2022 *SciPost Phys.* **13** 117

[52] Tang W and Khaymovich I M 2022 *Quantum* **6** 733

[53] Tarzia M 2022 *Phys. Rev. B* **105** 174201

[54] Skvortsov M A, Amini M and Kravtsov V E 2022 *Phys. Rev. B* **106** 054208

[55] Cai X, Lang L J, Chen S and Wang Y 2013 *Phys. Rev. Lett.* **110** 176403

[56] DeGottardi W, Sen D and Vishveshwara S 2013 *Phys. Rev. Lett.* **110** 146404

[57] Liu F, Ghosh S and Chong Y D 2015 *Phys. Rev. B* **91** 014108

[58] Das A K and Ghosh A 2022 *Phys. Rev. E* **105** 054121

[59] Ahmed A, Ramachandran A, Khaymovich I M and Sharma A 2022 *Phys. Rev. B* **106** 205119

[60] Lee S, Andreanov A and Flach S 2023 *Phys. Rev. B* **107** 014204

[61] Venturelli D, Cugliandolo L F, Schehr G and Tarzia M 2023 *SciPost Phys.* **14** 110

[62] Das A K and Ghosh A 2023 *J. Phys. A: Math. Theor.* **56** 495003

[63] Buijsman W 2024 *Phys. Rev. B* **109** 024205

[64] Čadež T, Kumar Nandy D, Rosa D, Andreanov A and Dietz B 2024 *New J. Phys.* **26** 083018

[65] Le Doussal P 2025 *J. Phys. A: Math. Theor.* **58** 335002

[66] Baron J W 2025 *Europhys. Lett.* **151** 21002

[67] Jahnke V, Nandy P, Pal K, Camargo H A and Kim K Y 2025 Free probability approach to spectral and operator statistics in Rosenzweig-Porter random matrix ensembles (arXiv:2506.04520)

[68] Lunkin A and Tikhonov K 2025 *SciPost Phys.* **19** 015

[69] Safonova E, Lunkin A and Feigel’man M 2025 Density of states correlations in Lévy Rosenzweig-Porter model via supersymmetry approach (arXiv:2505.23903)

[70] Rosenzweig N and Porter C E 1960 *Phys. Rev.* **120** 1698–714

[71] Sarkar M, Ghosh R, Sen A and Sengupta K 2021 *Phys. Rev. B* **103** 184309

[72] Sarkar M, Ghosh R and Khaymovich I M 2023 *Phys. Rev. B* **108** L060203

[73] Ghosh R, Sarkar M and Khaymovich I M 2025 *Phys. Rev. B* **111** L220102

[74] Wishart J 1928 *Biometrika* **20A** 32

[75] Bun J, Bouchaud J P and Potters M 2017 *Phys. Rep.* **666** 1–109

[76] Bouchaud J P 2017 Random matrix theory and (big) data analysis *Stochastic Processes and Random Matrices: Lecture Notes of the Les Houches Summer School* vol 104 (Oxford University Press)

[77] Semerjian G 2024 *J. Stat. Phys.* **191** 139

[78] Ledoit O and Péché S 2010 *Probab. Theory Relat. Fields* **151** 233–64

[79] Bun J, Allez R, Bouchaud J P and Potters M 2016 *IEEE Trans. Inf. Theory* **62** 7475–90

[80] Biroli G and Mézard M 2023 *J. Stat. Mech.* **093402**

[81] Marčenko V A and Pastur L A 1967 *Math. USSR Sb.* **1** 457–83

[82] Voiculescu D V, Dykema K J and Nica A 1992 *Free Random Variables* (American Mathematical Society) (available at: https://books.google.it/books?id=hzjd1KxaRroC)

[83] Zee A 1996 *Nucl. Phys. B* **474** 726–44

[84] Mirlin A D 2000 *Phys. Rep.* **326** 259–382
[85] Mott N and Twose W 1961 *Adv. Phys.* **10** 107

[86] De Luca A, Scardicchio A, Kravtsov V E and Altshuler B L 2013 Support set of random wave-functions on the Bethe lattice (arXiv:1401.0019)

[87] García-Mata I, Martin J, Dubertrand R, Giraud O, Georgeot B and Lemarié G 2020 *Phys. Rev. Res.* **2** 012020

[88] Gradshteyn I S and Ryzhik I M 2007 *Table of Integrals, Series and Products* 7th edn (Elsevier/Academic)

[89] David H Asymptotic expansion of integral involving an ArcTan Mathematics Stack Exchange (available at: https://math stackexchange.com/q/4523130)

[90] Mézard M, Parisi G and Virasoro M 1987 *Spin Glass Theory and Beyond (Lecture Notes in Physics Series) (World Scientific)* (available at: https://books.google.it/books?id=ZIF9QgAACAAJ)

[91] Edwards S F and Jones R C 1976 *J. Phys. A: Math. Gen.* **9** 1595-603

[92] Metz F L and Pérez Castillo I 2016 *Phys. Rev. Lett.* **117** 104101

[93] Metz F L 2017 *J. Phys. A: Math. Theor.* **50** 495002

[94] Pérez Castillo I and Metz F L 2018 *Phys. Rev. E* **97** 032124

[95] Vivo P 2020 *J. Phys. A: Math. Theor.* **54** 025002

[96] Dyson F J 1962 *J. Math. Phys.* **3** 140-56

[97] Pandey A 1995 *Chaos Solitons Fractals* **5** 1275-85

[98] Bru M F 1991 *J. Theor. Probab.* **4** 725-51

[99] Zavatone-Veth J A and Pehlevan C 2023 *SciPost Phys. Core* **6** 026

[100] Budnick B, Forer P, Vivo P, Aufiero S, Bartolucci S and Caccioli F 2025 *J. Phys. A: Math. Theor.* **58** 325001