Similarity between a many-body quantum avalanche model and the ultrametric random matrix model

In the field of ergodicity-breaking phases, it has been recognized that quantum avalanches can destabilize many-body localization at a wide range of disorder strengths. This has in particular been demonstrated by the numerical study of a toy model, sometimes simply called the ‘’avalanche model” or the ‘’quantum sun model” [Phys. Rev. Lett. 129 , 060602 (2022)], which consists of an ergodic seed coupled to a perfectly localized material. In this paper, we connect this toy model to a well-studied model in random matrix theory, the ultrametric ensemble. We conjecture that the models share the following features. 1) The location of the critical point is predicted sharply by analytics. 2) On the localized site, both models exhibit Fock space localization. 3) There is a manifold of critical points. On the critical manifold, the eigenvectors exhibit nontrivial multifractal behaviour that can be tuned by moving on the manifold. 4) The spectral statistics is intermediate between Poisson statistics and random matrix statistics, also tunable on the critical manifold. We confirm numerically these properties.

Of equal importance is to understand the boundaries of quantum thermalization and the conditions under which ergodicity breaking phase transitions may take place.At the current stage of research, it is particularly important to establish toy models of ergodicity breaking phase transitions, which exhibit clear features of critical behavior already in finite systems amenable to numerical simulations.Advances in this subject may bring new impetus for ongoing experimental activities [17][18][19][20][21][22][23][24][25], as well as provide new perspective into many-body localization [26][27][28][29], for which different perspectives about its stability in the thermodynamic limit have recently been formulated .
Avalanche theory has been studied mostly within toy models, one of which was called the "Quantum Sun Model" (QSM) in [70,72], which we will mostly refer to below.At this moment, models like the QSM are the only examples of a many-body ergodicity breaking transition in Hamiltonian systems that have the potential to be truly well understood, as they seem well accessible both analytically and numerically.In this paper, we want indeed to start a more detailed investigation of the QSM.The main Leitmotiv we use here, is the comparison of the QSM to a well-known model in random matrix theory, namely the "ultrametric model" (UM) [75].Despite the fact that the latter model is not commonly thought of as a many-body model, we conjecture that the behaviour of the QSM and the UM is, when cast in the same language, nearly identical in the thermodynamic limit for a broad range of model parameters.In particular, we expect the following similarities, also sketched in Fig. 1.
1.Both models have a natural parameter that we call α, and that the transition occurs at α c = 1/ √ 2.

FIG. 1.
Phase diagram of the avalanche models studied in this work, as a function of the tuning parameter α that drives the system from an ergodic phase (in which short-range spectral statistics comply with random matrix theory predictions) to a nonergodic phase (that exhibits Fock space localization).At the critical point, α = αc, the spectral statistics may not comply with random matrix theory predictions and the Hamiltonian eigenstates exhibit multifractality in the basis of uncoupled spin-1/2 particles.
These similarities would be in particular remarkable because the QSM, being a many-body model, has a number of disorder variables that scales like the physical volume of the system, hence like the logarithm of the Hilbert space dimension.In contrast, the UM has a number of disorder variables that scales like the Hilbert space dimension itself.Yet, in a certain sense, the UM is a more natural model than the QSM to describe the coupling of a small ergodic region to a perfectly localized system; whereas the QSM artificially pretends that the coupling is only to single spins S i of the localized system, the UM describes the more realistic [76] case where this particular coupling also involves all spins S j , j < i, closer to the ergodic region.
In this paper, we proceed concretely as follows.In Sec.II we introduce the two models under consideration, the QSM and the UM.We then introduce the indicators of the critical point in Sec.III, i.e., the level spacing ratio, participation and entanglement entropies, and the Schmidt gap.We establish two criteria of the critical point: (a) the level spacing ratio and the entanglement entropy of the most weakly coupled spin exhibit a scale invariant point, and (b) the first derivatives of the participation and entanglement entropies exhibit a sharp peak.In both models, these scale invariant points and the peaks of the derivatives almost perfectly coincide with the analytically predicted value for the critical point.In Sec.IV we then study Fock space localization on the nonergodic side and multifractal properties at the manifold of critical points.To this end, we extract the fractal dimension from the scaling of participation entropies.We argue that the degree of multifractality can be tuned by either modifying the size of the initial ergodic seed, or by the overall coupling of the ergodic seed to the remainder of the system.Finally, in Sec.V we study properties on the critical manifold through the lens of level statistics.We show that the values of the scale invariant level spacing ratio is close to the random matrix theory prediction if the fractal dimension is close to 1, and it decreases towards to Poisson value upon enhancing multifractality, i.e., upon decreasing the fractal dimension towards zero.We conclude in Sec.VI.

II. MODELS
In this section, we introduce the two models studied in this work, the QSM and the UM.In particular, we provide the definition of both models, norms of operators contained in the models, and in Appendix A we show their density of states.
In both models N denotes the number of spin-1/2 particles within the subsystem that we denote the quantum dot.Throughout our scaling analysis, we keep N = O(1) constant.Setting N = 3 as an example, the quantum dot is sketched as a region bounded by a blue circle in Fig. 2.However, the type of interaction that gives rise to the proceeding avalanche, as sketched in Figs.2(b)-2(d), depends on the particular model.With the term avalanche we have in mind a process that gives rise to ergodicity of spin-1/2 particles outside the dot.In the absence of coupling to the particles within the dot, the particles outside the dot are frozen in either of the S z j = ±1/2 states.The number of particles outside the dot is denoted by L and we consider the thermodynamic limit of both models by fixing N and sending L → ∞.

A. Quantum sun model (QSM)
The QSM Hamiltonian is given by The first term on the r.h.s. of Eq. ( 1), H dot , acts nontrivially on the dot degrees of freedom only.Properties of the latter are described by a 2 N × 2 N random matrix One of the main messages of this work is our conjecture that the UM of the random matrix theory (RMT) may be considered as a toy model of the avalanche theory, which retains all of the relevant physical features of the QSM while striping down all the unnecessary bits.Several aspects of the UM have been studied in the past [75,[77][78][79][80][81][82][83][84], however, its connection to the avalanche theory has to our knowledge not yet been explored.
We note that the initial motivation for studying the UM was to better understand the Anderson localization transition of noninteracting particles [75].Consequently, the term "ultrametric" corresponds to the particular geometry of a hierarchical lattice ("ultrametric lattice"), in which the metric is defined as the number of steps needed for a single particle to hop from one node to another node.Here we are interested in many-body physics and hence we formulate the UM in the Fock space of N +L spin-1/2 particles with dimension D = 2 N +L ; see the definition below.For simplicity we use the same name for the model.One may define an ultrametric distance d in our manybody model by ordering the spin-1/2 particles, cf.Fig. 2, by a decreasing coupling from the ergodic quantum dot.For simplicity, let us assume that there are no particles within the dot, i.e., N = 0.Then, the two spin configurations |σ The model Hamiltonian is constructed as a sum of block-diagonal random matrices Ĥk with k = 0, 1, . . ., L, where k is proportional to the ultrametric distance d.At each k, the matrix structure of Ĥk consists of 2 L−k diagonal blocks of size 2 N +k × 2 N +k .In analogy with the QSM, we sample each random block independently from the GOE distribution and use normalization in accordance with Eq. ( 2), such that Here, the superscript i denotes the i-th random block.We sample its matrix elements in analogy to the QSM, hence The full Hamiltonian of the UM then reads The first term Ĥ0 of size 2 N × 2 N models the initial quantum dot that consists of N spin-1/2 particles in the absence of coupling to the external particles.The sum in the second term mimics the exponentially decaying coupling between the dot and the k-th localized spin-1/2 particle through the exponentially decaying values of α k .Additionally, we have also included the parameter J tuning the overall perturbation strength, which carries certain analogies with the parameter g 0 in the QSM in Eq. (1).As in the QSM, we interpret the UM as being defined in a Fock space of N + L spin-1/2 particles, with D = 2 N +L .The matrix structure of the UM Hamiltonian at N = L = 3, α = 0.85 and J = 1 is shown in Fig. 3(b), and the density of states in both models are shown in Appendix A.
The UM can also be thought of as a discretized version of the power-law random banded model (PLRBM) [83,85,86].In the latter, the off-diagonal matrix elements h i,j are random numbers whose standard deviation decays with the distance r from the diagonal, at large r, as std(h |i−j|=r ) ∝ r −a = 2 −a log 2 r .The connection with the UM can be made if one approximates log 2 r with an integer, e.g., as d = ⌈log 2 r⌉, and associates d with the ultrametric distance of the UM.Then, the decay of the off-diagonal matrix elements in the PLRBM, ∝ 2 −ad , should be compared with the decay in the UM at large d, ∝ α d / √ 2 d , cf.Eqs. ( 3) and ( 4).This yields the relationship between the decay exponent a of the PLRBM and α in the UM, It follows from Eq. ( 5) that the transition point a c = 1 in the PLRBM corresponds to the transition point α c = 1/ √ 2 in the UM.A common feature of both models is the existence of 2 L = 8 dense blocks of size 2 N × 2 N = 8 × 8, denoting the ergodic quantum dot.Outside these blocks, the matrix is sparse in (a) and dense in (b).The matrix elements in each dense block in (a) are identical, except for their diagonal values, which are related to the distributions of random fields hj, see Eq. ( 1).The matrix elements of each block in (b) are sampled independently.When the size of the blocks is increased, the magnitude of their matrix elements is reduced accordingly.

C. Relationship between the QSM and UM
The matrix elements of both models, shown in Fig. 3 at N = L = 3, illustrate certain similarities and differences between the two models.The matrix elements within the dot are in both cases depicted by the dense blocks of size 2 N × 2 N = 8 × 8 along the diagonals.Then, the QSM only contains two-body interactions between particles outside the dot and randomly chosen particles within the dot.This gives rise to the overall sparse structure of the Hamiltonian matrix, shown in Fig. 3(a).In the UM, on the other hand, a particle k outside the dot (k = 1, 2, ..., L) is effectively coupled to all particles inside the dot as well as to particles k ′ < k outside the dot.This gives rise to the dense Hamiltonian matrix, shown in Fig. 3(b).
The most important similarity of both models is that the critical point α c between an ergodic and nonergodic phase is expected to occur in the thermodynamic limit (N fixed and L → ∞) at the same value, This expectation is based on the hybridization condition argument [62], which considers the situation when a particle k outside the dot interacts via a coupling g k with the ergodic bubble consisting of N + (k − 1) particles.
The hybridization condition is expressed as where the states |n⟩, |m⟩ are tensor product states of eigenstates of the ergodic bubble and the two-level system of the spin-1/2 particle k.Conjecturing that the matrix elements of eigenstates of the ergodic bubble satisfy eigenstate thermalization hypothesis [1,[5][6][7], one can approximate the matrix element as ⟨n| V |m⟩ ≈ g k √ ∆, where the level spacing ∆ of the many-body spectrum scales as ∆ ∝ 2 −(N +k) .In both models, the coupling g k scales as g k ∝ α k , such that the key tuning parameter is α ∈ [0, 1).The critical point then emerges when G ∝ α k 2 k/2 remains a nonzero constant when k → ∞, which occurs at α = α c = 1/ √ 2, as given by Eq. ( 6).We will show in this work that for system sizes amenable to state-of-the-art exact diagonalization of Hamiltonian matrices, Eq. ( 6) provides an accurate estimate of the transition point for a wide range of model parameters.In Sec.IV we will also consider an example of weak coupling between the inner and outer particles in the QSM, for which the transition point in the system sizes under investigation occurs at α that exceeds α c from Eq. (6).

III. INDICATORS OF THE CRITICAL POINT
We now turn our attention to the signatures of the critical point upon tuning the parameter α.We numerically study the statistical properties of the energy spectra and the properties of their corresponding energy eigenstates using full exact diagonalization of the Hamiltonian matrices from Eqs. (1) and (4).In all cases under investigation, we will first focus on the UM, for which finite-size convergence of the numerical results is the most favorable, followed by the comparison with the QSM.
The goal of this section is two fold.First, we would like to establish similarities between the two models under investigation, and show evidence that the prediction for the critical point from Eq. ( 6) provides an accurate estimate for the critical point in a broad range of parameters of both models.Then, we will compare different indicators of the critical point to explore which of them provide the most accurate prediction of the critical point, and to study for which regimes of model parameters the finite-size effects around the critical point scale most favorably.The latter will be instrumental for the analysis of Fock space localization in the nonergodic phase and multifractal properties at the critical point, which we will study in Sec.IV.

A. The average ratio of the adjacent level spacings
We first proceed by analysing the average ratio of the adjacent level spacings r [28], shortly the average gap ratio.For a target eigenlevel n, the ratio rn is defined as where δE n = E n+1 − E n is the level spacing between the eigenlevels n + 1 and n, respectively, while r n is the ratio of the consecutive level spacings, By definition, {r n } only assume values from the interval [0, 1], and hence no unfolding procedure is needed to eliminate the influence of finite-size effects through the local density of states.To obtain the average value r, we first average over N eig = 500 eigenstates near the center of the spectrum for each Hamiltonian realization, denoted by ⟨• • • ⟩ n in Eq. ( 9), and then over an ensemble of spectra for different Hamiltonian realizations, denoted by ⟨• • • ⟩ H .In the ergodic regime, r assumes the GOE value r GOE ≈ 0.5307 [87], while the prediction for energy levels with Poissonian statistics is r Poisson = 2 ln 2 − 1 ≈ 0.3863 [28].
The results for r versus α are shown for the UM (at J = 1) in Fig. 4 and for the QSM (at g 0 = γ = 1) in Fig. 5.There are two main messages from these results.
The first result is that the critical point is rather accurately predicted for both models by α c = 1/ √ 2 from Eq. (6).The location of the critical point is estimated by a crossing point or r versus α at different values of L, which signals scale invariance of r at the critical point.We will show in the next sections that this approach to detect the critical point agrees very well with those based on the analysis of participation and entanglement entropies.
The second result is that by increasing N the spectral statistics at the critical point gets closer to the GOE value r GOE .This is most prominent at N = 3 in the UM, see Fig. 4(c), and at N = 5 in the QSM, see Fig. 5(b), for which the point α = α c signals the breakdown of r from r GOE .In these cases, therefore, the criterion for the critical point based on the breakdown of r from r GOE replaces the criterion based on the emergence of a scale invariant crossing point of the r values.
Figure 5 also reveals an important numerical detail about the QSM.Namely, at N = 3 the crossing point of the r values emerges at α that is slightly larger than α c , see Fig. 5(a), while at N = 5 the breakdown point of r from r GOE is accurately predicted by α c .This property was already noticed in Ref. [70] that considered N = 3 and in Ref. [73] that considered N = 5.We conclude from these results that by increasing N in the QSM, the transition point for the available system sizes gets more accurately predicted by α c from Eq. (6).Still, since the thermodynamic limit is obtained by increasing L → ∞ and hence L is expected to be larger than N , the value of N in the actual calculations should not be too large.In what follows, we set N = 5 as a compromise to obtain optimal convergence of finite-size effects in the QSM.
In Figs.16 and 17 of Sec.IV we will further show that decreasing J and g 0 in the UM and the QSM, respectively, further shifts the r values at the critical point closer to the Poisson value r Poisson .This insight will be used to tune the multifractal properties of the eigenstates at the critical point.

B. The inverse participation ratio and the participation entropy
We next focus our attention towards the properties of the Hamiltonian eigenstate wavefunctions.To that end, we calculate the inverse participation ratio (IPR) for a selected number of Hamiltonian eigenstates |n⟩ .The IPR, generalized to an arbitrary index q, is defined as The generalized IPR as such is basis dependent.Here, |i⟩ are the states in the computational basis, i.e, the basis in which Ŝz i operators are diagonal, |n⟩ is a Hamiltonian eigenstate at the corresponding energy E n , and q is a parameter of the calculation.
For convenience, the quantity that we actually study is the participation entropy, which is defined as the logarithm of the IPR.Here we mainly focus on the participation entropy of the typical IPR, where the brackets ⟨• • • ⟩ n and ⟨• • • ⟩ H denote the averages over Hamiltonian eigenstates at a fixed realization and different Hamiltonian realizations, respectively.The limiting value of S (q) typ as q → 1 is the von Neumann participation entropy, We note that we have also studied the participation entropy of the average IPR, S (q) = 1 1−q ln⟨⟨P −1 q ⟩ n ⟩ H , with no qualitative physical differences (not shown here).The goal of this section is to explore to which extent one can use the IPR-based quantities to pinpoint the critical point.Figures 6(a)-6(c), and the corresponding insets, show S (q) typ versus α at various q in the UM at J = 1 and N = 1.They behave according to the expectations: S (q) typ approaches in the limit α → 0 an L-independent constant, and it increases (roughly linearly) with L in the ergodic phase at α > α c .However, based on the analysis of S (q) typ only, it appears to be impossible to determine the critical point without any a priori knowledge of it.
In Figs.6(d)-6(f) we show that a valuable information about the critical point can be obtained by calculating the derivative of the participation entropy, dS (q) typ /dα, and plot it versus α.This analysis is inspired by the recent results for the Anderson models [88,89], which showed that the peak of dS (q) typ /dW (where W denotes the disorder amplitude) is located very close to the Anderson localization transition in three dimensions [88], and it drifts to zero in two dimensions [89].A similar tendency is also observed in the UM at J = 1 and N = 1, shown in Figs.6(d)-6(f).Namely, the peak in dS (q) typ /dα is, for most values of q, located very close to the predicted critical point α = α c from Eq. ( 6).
As a technical remark, we calculate the derivatives dS (q) typ /dα by first interpolating the raw data for S (q) typ with cubic splines using the UnivariateSpline function from the scipy.interpolatelibrary.Then, using the tools available within the same interpolating function, we also evaluate the first derivatives.We use the same procedure also to evaluate the derivatives of the Rényi entanglement entropies in Sec.III C.
The q-dependence of the results in Figs.6(d)-6(f) reveals rich information about the behavior of participation entropies close to the critical point.We observe typ [panels (a-c)] and their derivatives dS (q) typ /dα [panels (d-f)] in the UM at J = 1 and N = 1.(a) and (d): q = 0.5 in the main panel and q = 0.1 in the insets.(b) and (e): q = 1.(c) and (f): q = 2 in the main panel and q = 6 in the inset.Vertical dashed lines denote α = αc from Eq. ( 6).Results for S (q) typ are averaged over N samples = 1000, 400, 300 Hamiltonian realizations for L ≤ 12 , L = 13, L = 14, respectively.Averaging over the eigenstates, ⟨• • • ⟩n, was performed over Neig = 1000 eigenstates near the center of the spectrum for each data point.
typ [panels (a-c)] and their derivatives dS (q) typ /dα [panels (d-f)] in the QSM at g0 = γ = 1 and N = 5.(a) and (d): q = 0.5 in the main panel and q = 0.1 in the insets.(b) and (e): q = 1.(c) and (f): q = 2 in the main panel and q = 6 in the inset.Vertical dashed lines denote α = αc from Eq. ( 6).Results for S (q) typ are averaged over N samples = 500 Hamiltonian realizations.Averaging over the eigenstates, ⟨• • • ⟩n, was performed over Neig = 500 eigenstates near the center of the spectrum for each data point.
that at q = 2, see Fig. 6(f), the peak location coincides with α = α c to extremely high accuracy even in small systems.This remarkable property is also observed at higher q, see the inset of Fig. 6(f), while at smaller q this high precision is lost, see the results in Figs.6(d)-6(e).
In the localized regime of the PLRBM (i.e., at a > 1), it was shown that the generalized IPR at q < 1/2 may not exhibit localization at arbitrary a > 1 as a consequence of power-law localization [85,86].Based on the relationship between the PLRBM and UM discussed in Sec.II B, we expect similar features to emerge in the UM as well.This is consistent with the observation in the inset of Fig. 6(d), which suggests that at q ≪ 1 in the UM it becomes nearly impossible to detect the critical point since dS typ /dα exhibits a shoulder in the vicinity of the critical point rather than a sharp peak.
While the results in Fig. 6 were obtained for the UM at N = 1, in Fig. 20 of Appendix B we also show the results (at q = 2) for other values of N .We observe that the finite-size effects scale most favorably at N = 1, i.e., at this value of N the peak of dS (q) typ /dα emerges almost exactly at α = α c already in small systems of N +L = 11 spin-1/2 particles.We hence consider the UM at N = 1 as the model that is best suited for the analysis of the critical behavior in finite systems, and we only study the N = 1 case further on.
In Fig. 7 we complement the results for the participation entropies in the UM by showing the results for the QSM at g 0 = γ = 1 and N = 5.Qualitatively, the results in Fig. 6 and 7 are rather similar, which is the main message of these analyses.Still, the accuracy to determine the location of the critical point via the peak of dS (q) typ /dα is in the QSM not as high as in the UM. Figure 7(f) shows that again the most accurate results are obtained in the large q regime, q ≳ 2. In case the precise position of the critical point is not known in advance, a possible method to improve the prediction for the critical point (that we do not pursue here) is to extract the position of the peak for each L, and then scale these values to the limit L → ∞.
To summarize the results for the participation entropies, we note that its derivative w.r.t. the tuning parameter of the transition, dS (q) typ /dα, exhibits a peak that we expect is located at the critical point in the thermodynamic limit.An additional insight that we obtained from our analysis is that in finite systems amenable to exact diagonalization, the most accurate prediction for the critical point are obtained at large values of q ≳ 2. On the other hand, the results in the opposite limit q ≪ 1 are, at least for the available system sizes, not very useful for determining the critical point.

C. Rényi entanglement entropies of eigenstates
While the analysis of the participation entropies in the previous section turned out to be quite useful for the determination of the critical point, we here complement these results by studying the entanglement based mea-sures.The participation entropies contain information about the whole wavefunction, and hence they may not be very sensitive to the properties of the particles that are "most distant" from the ergodic quantum dot.(With "most distant" we have in mind the particles that are most weakly coupled to the ergodic quantum dot.)On the other hand, the entanglement entropies allow for selecting arbitrary subsystems.Here we focus on the entanglement properties of most distant spin-1/2 particles.A physical motivation for this choice is that the ergodicity breaking phase transition within the avalanche theory is expected to occur when the avalanche fails to thermalize the most distant particles, and hence the change of their properties contains crucial information about the transition.
We consider the eigenstate entanglement entropies of subsystems that consist of p most distant spin-1/2 particles.The Fock space H of the studied Hamiltonian carries a tensor product structure where H N refers to the dot degrees of freedom and H i to the i-th spin-1/2 particle outside the dot.We partition the system into two partitions A p and B p , such that H = H Ap ⊗H Bp .Here, we take B p to be a collection of p spins farthermost from the dot, while A p denotes the remainder of the system that consists of N + L − p particles.In the following, unless necessary, we shall omit the subscript p when referring to the partitions.For the density matrix ρ of the full system, we obtain the reduced density matrix of the subsystem B by tracing out the degrees of freedom in A, where ρ = |n⟩⟨n| is the density matrix associated to the Hamiltonian eigenstate |n⟩.Upon diagonalization of ρB we obtain the eigenvalue spectrum of the reduced density matrix, which we denote by {λ 1 , . . ., λ i , . . ., λ D B }, where λ i ≥ λ i+1 .Here, D B is the reduced Fock space dimension which equals 2 p in our analysis.The q-th Rényi eigenstate entanglement entropy, which is normalized such that its maximum value is bounded from above by 1, is then defined as with q > 0 and q ̸ = 1.The limiting value of S (q) as q → 1 is the von Neumann eigenstate entanglement entropy, 0.0 0.5 1.0 (a) q = 0.5 0.5 0.9  q) [panels (a-c)] and their derivatives dS (q) /dα [panels (d-f)] in the UM at J = 1, N = 1 and p = 1.(a) and (d): q = 0.5 in the main panel and q = 0.1 in the insets.(b) and (e): q = 1.(c) and (f): q = 2 in the main panel and q = 6 in the inset.Vertical dashed lines denote α = αc from Eq. ( 6).Results for S (q) are averaged over N samples = 1000, 400, 300 Hamiltonian realizations for L ≤ 12 , L = 13, L = 14, respectively.Averaging over the eigenstates, ⟨• • • ⟩n, was performed over Neig = 1000 eigenstates near the center of the spectrum for each data point.q) [panels (a-c)] and their derivatives dS (q) /dα [panels (d-f)] in the QSM at g0 = γ = 1, N = 5 and p = 1.(a) and (d): q = 0.5 in the main panel and q = 0.1 in the insets.(b) and (e): q = 1.(c) and (f): q = 2 in the main panel and q = 6 in the inset.Vertical dashed lines denote α = αc from Eq. ( 6).Results for S (q) are averaged over N samples = 500 Hamiltonian realizations.Averaging over the eigenstates, ⟨• • • ⟩n, was performed over Neig = 500 eigenstates near the center of the spectrum for each data point.
As in Eqs. ( 9), ( 11) and ( 12), the averaging is performed over Hamiltonian eigenstates and different Hamiltonian realizations.At p = 1 one can get further analytical insight into the entanglement entropies, which we present in Appendix C. Results for the eigenstate entanglement entropies S (q) at p = 1 are shown for the UM in Fig. 8.They are consistent with the limiting behaviors described in Eqs.(C2) and (C3) of Appendix C, i.e., S (q→0) is very close to 1 at α = α c , see the inset of Fig. 8(a), and S (q→∞) at α = α c is considerably smaller than 1, see the inset of Fig. 8(c).
The eigenstate entanglement entropies S (q) exhibit an important advantage when compared to the participation entropies S (q) typ , namely, the critical point can readily be estimated to high accuracy from the crossing point of S (q) versus α at different L, see Figs. 8(a)-8(c).Moreover, its derivative, dS (q) /dα, also exhibits a peak that is located very close to the predicted critical point α = α c .We again observe that the location of the peak almost exactly coincides with α c at large q, see Fig. 8(f), while at small q the agreement is less accurate, see Fig. 8(d).To evaluate the derivatives of the entanglement entropies, we use the same procedure as the one outlined in Sec.III B.
Analogous results for the eigenstate entanglement entropies S (q) at p = 1 are shown for the QSM in Fig. 9. Also in this case, there exist a scale invariant (i.e., L independent) point of S (q) that is almost exactly located at the critical point α = α c , see Figs. 9(a)-9(c).The derivative dS (q) /dα exhibits a peak very close to the critical point, see Figs. 9(d)-9(f).Nevertheless, the location of the peak in the QSM does not coincide with α = α c as accurately as in the UM shown in Figs.8(d)-8(f).Still, the common feature of both models is that the position of the peak gets closer to α c when q is increased.All these results confirm that the similarity of the critical behavior in both models also emerges on a level of entanglement entropies.
So far, we mostly focused on the q dependence of S (q) and its impact on the determination of the critical point.Results in Figs. 8 and 9 showed that one can pinpoint very accurately the critical point using S (q) at p = 1, i.e., when studying the entanglement properties of the subsystem that consists of a single, most distant spin-1/2 particle from the ergodic quantum dot.A natural question is then to ask what is the optimal size of the subsystem that is most sensitive to the emergence of the critical behavior.In Figs.21 and 22 of Appendix C we show results that are analogous to those in Figs. 8 and 9, respectively.However, they are calculated for the subsystems that consist of p = 4 most distant particles from the dot.Results in Figs.21 and 22 suggest that they are still consistent with the emergence of the critical point in the vicinity of α = α c , nevertheless, the accuracy is not as good as for the results at p = 1, in particular for the QSM.From these we conclude that, at least for the QSM, the most valuable information about the transition is encoded in the particles that are most distant from the ergodic quantum dot.In other words, when the ergodic bubble fails to thermalize the entire system, these particles are the first to exhibit nonergodic properties.

D. Schmidt gaps
As another indicator of the critical point, which is related to the entanglement entropies, we consider the Schmidt gap ∆.The latter is defined as the difference between the two largest eigenvalues of the reduced density matrices of the subsystem B from Eq. ( 14), ∆ = λ 1 − λ 2 .In the actual numerical calculations, we then average ∆ over the midspectrum eigenstates and over different Hamiltonian realizations, similar to the other quantities studied before.Also similar to the other quantities, in ∆ we omit the index p indicating the number of spin-1/2 particles in the subsystem.
As already discussed in the previous section, it is convenient to study the case of a single two-level system (p = 1), for which λ 2 = 1 − λ 1 and hence Results for ∆ at p = 1 are shown for the UM and the QSM in the main panels of Figs.10(a) and 11(a), respectively.They provide an extremely accurate measure of the critical point.In particular, we observe that ∆ is scale invariant at α = α c , and it tends towards ∆ → 0 at α > α c and ∆ → 1 at α < α c in the thermodynamic limit L → ∞.This establishes the Schmidt gap ∆ as a candidate that may play a role of an order parameter for the ergodicity breaking phase transition.A natural question is to ask about the optimal size of the subsystem that exhibits most sharp signatures of the transition.To this end we extend the analysis of ∆ to p = 4 in Figs.10(b) and 11(b), i.e., to the case where the subsystem consists of four most distant particles from the ergodic quantum dot.Also in this case, one can still detect the critical point by inspecting the position of the scale invariant point of ∆.However, in this case the scale invariant value of ∆ is much closer to zero, and hence the signatures of the scale invariance are not as sharp as in the case of p = 1.
At p = 1, both the Schmidt gap ∆ and the entanglement entropies S (q) are functions of a single eigenvalue λ 1 , as expressed by Eqs. ( 18) and (C1), respectively.It is then expected that the S (q) is a well-defined function of ∆, independent of the system size L.This property is demonstrated in the insets of Figs.10(a) and 11(a) at q = 2.An interesting detail of these figures, though, is that the results do not exactly follow the expected behavior that would have been derived from Eqs. ( 18) and (C1) if there was no averaging over eigenstates and Hamiltonian realizations.It then appears to be not entirely trivial that at p = 1, S (2) is a well-defined function of ∆ despite this function being different from Eq. ( 19).As a side remark, in Fig. 23(a) of Appendix D we show that when no averaging over Hamiltonian eigenstates and Hamiltonian realizations are performed, the second Rényi entanglement entropy is indeed a well-defined function of the Schmidt gap as predicted by Eq. ( 19).The deviations between Eq. ( 19) and the numerical results in the insets of Figs.10(a) and 11(a) are hence a consequence of the averaging.
At p > 1, we are not aware of any formal argument to predict a unique relationship between S (q) and ∆.Quite surprisingly, however, we still observe a nearly perfect collapse of the results for S (2) when plotted at p = 4 as a function of ∆, see the insets of Figs.10(b) and 11(b).A detailed inspection in Fig. 23 of Appendix D reveals that the collapse is not present when S (2) is plotted versus ∆ for a single Hamiltonian eigenstate, and that the signatures of a well-defined functional dependence of S (2) on ∆ can be readily observed after the averaging over eigenstates within a single Hamiltonian realization.While we are here not able to fully rationalize this behavior, we note that the remarkable scaling collapses share close similarities in both models.

IV. MULTIFRACTALITY AND FOCK SPACE LOCALIZATION
Having established in Sec.III the sharp agreement between the numerically determined location of the critical point and the analytical prediction, we here focus on the properties on the nonergodic side and at the critical point.The central quantity of study will be the fractal dimension introduced below.

A. Fractal dimension
In this work we extract the fractal dimension from the participation entropies, which were introduced in Sec.III B. In particular, in Sec.III B we focused on the participation entropies S (q) typ that we calculated from the typical IPR, see Eq. (11), and hence we refer to the corresponding fractal dimensions as d (q) typ .We calculate d (q) typ from the ansatz [90] where D = 2 N +L is the Fock space dimension and b (q) typ is an L-independent constant.The ansatz in Eq. ( 20) is phenomenological, and it is expected to describe well the results in the asymptotic regime, i.e., at sufficiently large L (the details of the numerical implementation will be discussed below).The fractal dimension d typ is an Lindependent constant that may depended on q.If that is the case, we refer to the wavefunction properties as being multifractal.We note that the fractal dimension is a basis-dependent quantity since it depends on the wavefunction coefficients c i = ⟨i|n⟩ calculated in some basis {|i⟩}.Here, {|i⟩} corresponds to the computational basis.
A similar but nonequivalent way to define the fractal dimension is via the average IPR, which we denote as ⟨⟨P −1 q ⟩ n ⟩ H .In this case, one extracts the decay coefficient τ (q) of the IPR from the ansatz from which one can define the fractal dimension d (q) via the relation d (q) = τ (q) /(q − 1).Nevertheless, studies of single-particle Anderson localization transition noted that the distribution of the IPRs at the critical point may be broad [91], and hence a more natural choice of consideration is d (q) typ instead of d (q) .In all numerical results reported here we only focus on d (q) typ .We observe (not shown here) no significant differences between d (q) and d (q) typ , however, a quantitative comparison between d (q) and d (q) typ is left for future work.For the sake of completeness, we also introduce the decay coefficient of the typical IPR, denoted as τ Both quantities d (q) typ and τ (q) typ will be discussed below.The behavior of the fractal dimension is well understood in the ergodic phase and in the nonergodic phase that exhibits Fock space localization.These considerations apply to both d (q) and d (q) typ [92], and for simplicity we here only discuss the former.In the ergodic phase that exhibits many-body quantum chaos [1], the wavefunction coefficients c i = ⟨i|n⟩ can be considered as normally distributed random variables with zero mean and variance 1/D, and hence ⟨⟨P −1 q ⟩ n ⟩ H scales with the Hilbert space dimension approximately as giving rise to the decay coefficient τ (q) = q − 1 and a qindependent fractal dimension d q = 1.Specifically, in the commonly studied case at q = 2, the IPR scales as ⟨⟨P −1 2 ⟩ n ⟩ H ∝ D −1 , which is in agreement with the exact result predicted by the GOE of the RMT [93], ⟨⟨P −1 2 ⟩ n ⟩ H = 3/D.In the opposite case of Fock space localization, one may think of only O(1) coefficients c i begin nonzero, and thus ⟨⟨P −1 q ⟩ n ⟩ H does not scale with the system size.This implies τ (q) = d (q) = 0 for all q.The only exception is q = 1, at which P −1 q = 1 due to normalization, and hence d (1) cannot distinguish between the ergodic and nonergodic phase.However, studying the typical fractal dimension d (q) typ that is the focus of our study, no such limitation occurs at q = 1.
In finite systems such as those studied here, one often obtains d (q) typ that is away from the limiting cases discussed above, i.e., 0 < d (q) typ < 1.This result may either be consistent with (multi)fractality of the system, or suggesting that the system is not yet in the asymptotic regime.In our numerical analyses, we first verify whether the system can be considered as being in the asymptotic regime, and hence if the ansatz from Eq. ( 20) can be applied.Figures 12, 13 and 24 show examples in which the system is close to the asymptotic regime: these occur deep in the ergodic regime at α c ≪ α ≲ 1 and at the critical point α = α c .In contrast, outside the asymptotic regime, d (q) typ should be replaced by a number that depends on L, i.e., d typ , which we calculate in Fig. 12.We numerically extract d (q,L) typ by computing the values of slope of S (q) typ between consecutive system sizes L − 1 and L, d Such procedure was recently used in studies of multifractal properties of wavefunctions on different types of Anderson graphs [94][95][96].From the flow of d (q,L) typ with L, one may conjecture about the fate of d (q) typ in the thermodynamic limit L → ∞.

B. Fock space localization on the nonergodic side
We first focus on the properties on the nonergodic side of the critical point, at α < α c .Previous studies in the UM have established, almost up to a rigorous level, that in the nonergodic phase the system exhibits Fock space localization [75,78].
Here we are particularly interested on the nonergodic side of the QSM.The latter can be divided in two regimes: the regime close to the critical point and the regime deep in the nonergodic phase.Figure 12 shows the L-dependent fractal dimension d (q,L) typ from Eq. ( 24), which suggest that at sufficiently large q ≳ 1, the first regime approximately belongs for the interval 0.5 ≲ α < α c .In contrast, at α ≲ 0.5 and q ≳ 1 the system exhibits clear signatures of Fock space localization since d (q) typ ≈ 0. These results are corroborated in Figs. 14 and 15 by the results at α = 0.2 in the UM and the QSM, respectively, which consistently show d (q) typ ≈ τ (q) typ ≈ 0. In fact, we even observe slightly negative values of d (q) typ , see the results at α = 0.2 in Figs.13(a) and 15.In this case, we expect d (q) typ → 0 in the thermodynamic limit.The limit q → 0, which is not the focus of this study, may represent an exception to these considerations.As discussed in Sec.II B, it is understood from the studies of the PLRBM [85,86] that the observation of Lindependent IPR (and consequently the participation entropies) may only emerge deep in the localized regime, i.e. for sufficiently large a, with the threshold value for a depending on q.Similar arguments likely apply to the UM, and it is beyond the scope of this paper to establish to what degree these arguments also apply to the QSM.
The most interesting aspect of Fig. 12 is the flow of d (q,L) typ upon increasing L. In the nonergodic regime at 0.5 ≲ α < α c , the results are consistent with a flow towards Fock space localization in the thermodynamic limit, i.e., d (q,L) typ → 0 when L → ∞.These results are also consistent with the recent analysis of d (q) at α = 0.6 and q = 2 in the QSM [73], which was interpreted as Fock-space localization ergodic 0.5 α Fock-space localization ergodic 0.5 0.8 α Fock-space localization ergodic FIG.12.The L-dependent fractal dimension d (q,L) typ from Eq. (24) in the QSM at g0 = γ = 1 and N = 5, at (a) q = 0.1, (b) q = 1, (c) q = 2 and (d) q = 6.The arrows denote the flows towards Fock space localization in the asymptotic regime.Vertical dashed lines denote α = αc from Eq. ( 6).(a) q = 0.1 q = 0.3 q = 0.5 q = 0.7 q = 0.9 q = 1.5 q = 2.0 q = 3.0  20) for Fock space dimensions of the total system in the range from D = 2 13 to D = 2 16 .They allow for the extraction of the fractal dimension d (q) typ and the coefficient τ (q) typ , plotted in Fig. 15.
On the other hand, the results at α > α c in Fig. 12 exhibit values that are close to, but not necessary equal to d (q,L) typ = 1.In the vicinity of the critical point, d (q,L) typ is nonzero and lower than 1, which we interpret as a signature of multifractality and will be discussed in more details in the next section.
Summarizing these results, we expect that Fock space localization is a property of the entire nonergodic phase in both the UM and the QSM.This is different from what is expected to occur in the putative MBL phase, for which emergence of multifractality was conjectured, based on simple theoretical arguments, in the entire phase [46,90,[97][98][99][100].

C. Multifractality on the critical manifold
We now focus our attention on the eigenstate wavefunction properties at the critical point.In the UM it is understood that the wavefunction is multifractal [77,80,83].Specifically, it was shown that d (q) exhibits q dependence and in the limit J → 0, explicit expressions were derived for d (q) .In lowest order in J, the q dependence of d (q) is the same for the UM and the PLRBM [77,85].
Beyond the UM, multifractality was observed in several random matrix ensembles [91,[101][102][103].Its concept was also extended to physical models without interactions such as Anderson models on hypercubic lat- typ versus q, and (b) the decay coefficient τ typ of the typical IPR q, in the UM at N = 1.The values are extracted from the S (q) typ versus ln D curves, such as those shown in Fig. 24, using the ansatz in Eq. ( 20).The dashed lines in both panels correspond to the ergodic limit d (q) typ = 1, τ (q) typ = q − 1, and the dashed-dotted lines correspond to the nonergodic limit d g 0 = γ = 1, α = 0.9 typ versus q, and (b) the decay coefficient τ (q) typ of the typical IPR versus q, in the QSM at N = 5.The values are extracted from the S (q) typ versus ln D curves, such as those shown in Fig. 13, using the ansatz in Eq. ( 20).The dashed lines in both panels correspond to the ergodic limit d (q) typ = 1, τ (q) typ = q − 1, and the dashed-dotted lines correspond to the nonergodic limit d tices [86,92,[104][105][106][107] and on graphs [94][95][96][108][109][110], and the connections were made to quantum dynamics [73,[111][112][113][114][115][116][117][118].A currently very active direction of research is to explore to which degree the concept of multifractality can be applied to many-body quantum systems at the boundary of quantum chaos, such as those studied in the context of MBL [46,90,[97][98][99][100].
Here we first complement previous results on multifractality in the UM by showing that also d (q) typ exhibits multifractal properties.Specifically, in Figs.24(b) and 24(c) in Appendix E we show that S (q) typ can be well fitted by the ansatz from Eq. ( 20), suggesting that the system is indeed very close to the asymptotic regime.We then show the extracted values of d (q) typ and τ (q) typ as a function of q in Fig. 14.The fractal dimension d (q) typ decreases with q, and this dependence on q confirms the multifractal character of the wavefunction.typ and τ (q) typ versus q for two different parameter values J = 1 and J = 0.5 at the critical point α = α c .Different values of J give rise to different fractal dimensions of wavefunctions, suggesting the emergence of a manifold of critical points.In Sec.V we further study the impact of J on the critical properties, specifically, on the spectral statistics.It is possible that there is a one-parameter family of critical points tuned by the parameter J, however, we cannot rule out the existence of a higher-dimensional critical manifold.
The central question is whether the QSM also exhibits multifractal properties at the critical point.If the answer is affirmative, the next question is then to explore whether one can tune the multifractal properties within the critical manifold as in the UM.
The results in Fig. 12 have already suggested that the L-dependent fractal dimension d (q,L) typ exhibits no, or only mild, dependence on L at the critical point of the QSM.This observation is confirmed by the results in Figs.13(b) and 13(c), which exhibit a linear dependence of S (q) typ on ln D according to Eq. ( 20), and hence suggest that the system is already very close to the asymptotic regime in which the L-dependence on d (q) typ should not be large.The results for the fractal dimension d (q) typ , and for the decay coefficient τ (q) typ of the typical IPR, are shown in Fig. 15.At the critical point, they exhibit two important features: typ is nonzero and lower than 1 for all nonzero values of q under consideration, and it exhibits clear q dependence (it monotonically decreases with q).While the first property may be argued as being consistent with the results of Ref. [73], the second property of the QSM has to our knowledge not yet been explored.We interpret these results as evidence of the multifractal character of the wavefunctions.
The results for the QSM in Fig. 15 appear to be very similar to those for the UM in Fig. 14.Interestingly, one also observes similarity when studying the role of the coupling parameters J and g 0 in the UM and the QSM, respectively.In both cases, tuning the parameters J and g 0 from large to small values gives rise to a decrease of d (q) typ , which we interpret as tuning the multifractality from weak to strong.
Note that the results in Sec.III have established that the location of the critical point in the QSM is sharply predicted by α = α c from Eq. ( 6) when all the model parameters are of the same order, which is certainly the case at g 0 = γ = 1.If, however, the coupling parameter g 0 is varied, which is going to be studied in more detail in Sec.V, the transition point observed in finite systems may not accurately coincide with α c .For example, at g 0 = 0.3 that is also studied in Fig. 15, the transition point from the r statistics appears to be located around α = 0.81, see Fig. 17(b).This is the reason why multifractal properties of the QSM in Fig. 15 are, at g 0 = 0.3, studied at α = α and not at α = α c .

V. SPECTRAL STATISTICS ON THE CRITICAL MANIFOLD
Having established multifractal properties at the critical point of both models, we now turn our attention to the spectral properties that have already been studied in Sec.III A. We are particularly interested in how the spectral properties at the critical point change when the system is tuned from weak multifractality, i.e., from large fractal dimension d (q) typ (at q ≈ 1) towards strong multifractality, i.e., towards small fractal dimension d (q) typ .In Figs.16 and 17 we study the behavior of the average gap ratio r, defined in Eq. ( 9), around the critical point at different values of model parameters.We focus on the value of r at the critical point, which is detected by the scale invariant crossing point of r versus α.While in the UM the crossing point occurs at the predicted value α = α c for essentially all model parameters of investigation, we already highlighted in the previous sections that this may not necessary be the case in the QSM.An example of the latter is given in Fig. 17(b), where the crossing point at g 0 = 0.3 emerges close to α = α = 0.81.Clarifying the fate of this crossing point in the thermodynamic limit appears to be a challenging task that is beyond the scope of the paper.We note that at this point we are not aware of any rigorous argument that would prevent the crossing point α from drifting towards the predicted critical point α c from Eq. ( 6).
Figure 16 studies the impact of the coupling parameter J in the UM on the nature of the critical point.From Fig. 14 we have already learned that decreasing J gives rise to smaller values of the fractal dimension d (q) typ and hence to stronger multifractality.Figure 16 suggests that stronger multifractality is associated with a smaller value of r at the critical point, i.e., with spectral statistics that is closer to the Poisson distribution.We hence expect that J represents the tuning parameter of the spectral statistics within the critical manifold, spanning from the RMT-like statistics at weak multifractality to Poissonlike statistics at strong multifractality.These properties share some analogies with the PLRBM [85], in which the characteristic length b tunes the properties of the critical point, ranging from Poisson-like at strong multifractality (b ≪ 1) to RMT-like at weak multifractality (b ≫ 1) [74,119].
Similar behaviour is observed in the QSM as a function of the coupling parameter g 0 , see Fig. 17.In particular, large g 0 > 1 drive the critical point towards RMT-like statistics, see Fig. 17(a), while small g 0 < 1 bring it closer to the Poission-like statistics, see Fig. 17(b).The latter is consistent with strong multifractality observed in Fig. 15 at g 0 = 0.3.On the other hand, the change of γ, which governs the spectral width of the ergodic quantum dot in the QSM, does not appear to have any significant impact on the spectral statistics within the critical manifold, see the insets of Fig. 17.

VI. CONCLUSION
This work establishes a direct connection of the critical behavior between the two different models, the UM of the RMT and the QSM.A convenient aspect of both models is that there exist quantitative analytical arguments for the value of the critical point of the ergodicity breaking phase transition.Carrying out exact numerical calcula- tions for various ergodicity indicators, we showed that these analytical arguments sharply predict the location of the critical point.While our main goal was to establish the similarity of the two models, we also introduced certain ergodicity measures that have previously not received much attention.Among those, we stress the entanglement entropy of the most distant (i.e., most weakly coupled) spin-1/2 particle, which exhibits scale invariant behavior at the critical point, and the derivatives of the participation and entanglement entropies, which exhibit a sharp peak at the critical point.Another feature of the numerical calculations is that the location of the critical point in the QSM is closest to the predicted analytical value when the critical properties comply with the RMT predictions, and it exhibits small deviations otherwise.While this feature was already noticed before [70,72,73], future work should explore in more details the fate of these small deviations in the thermodynamic limit.
On the nonergodic side of the transition, we argued that both models exhibit Fock space localization in the eigenbasis of the Ŝz operator.While the latter was already established in the UM [75,78], our results suggest that also in the QSM, Fock space localization may emerge in the entire nonergodic phase in the thermodynamic limit.This result is different from what was proposed for the putative MBL phase in random-field spin-1/2 chains [46,90,[97][98][99][100], in which the entire nonergodic phase was conjectured to be multifractal.Therefore, while one can draw certain parallels between the ergodicity breaking in the QSM and MBL, there are also clear differences between these two phenomena.
At criticality, the eigenvectors exhibit multifractal behaviour in the eigenbasis of the Ŝz operator, characterized by a family of nonzero fractal dimensions that are lower than unity.The fractal dimensions may vary with the size of the initial ergodic seed and the overall coupling of the ergodic seed to the remainder of the sys-tem, thereby suggesting the emergence of a manifold of critical points.This may carry some similarities with certain RMT-based models for Anderson localization transition, such as the power-law random banded matrix models [85] that exhibit a one-parameter family of critical points.However, whether there is a one-parameter family, or eventually a higher manifold of critical points in the QSM, is an open question that should be addressed in the future.
As a consequence of the emergence of the manifold of critical points, the spectral statistics on the critical manifold may be continuously varied from the RMT-like statistics to the Poisson-like statistics.In the former, the fractal dimensions are close to unity, while in the latter they are close to zero.In the context of singleparticle transitions, such as those in the power-law random banded matrices, this relationship is in certain limits understood even analytically [119], while a systematic study in many-body quantum systems appears to be a natural next step of investigation.
To summarize, our results reinforce the QSM and the UM as fertile playgrounds to study many-body ergodicity breaking phase transitions, and call for further characterization of their properties both from the numerical as well as from the analytical side.In Sec.II we introduced both models under investigation, the QSM and the UM.Here we provide additional information about the models by studying their coarsegrained density of states ρ(E) = δN/δE, which counts the number of Hamiltonian eigenstates δN in a narrow energy window of width δE.
Figure 18 shows the density of states ρ(E) in the QSM.It is accurately described by the normal distribution when g 0 = γ = 1, see the inset of Fig. 18, while deviations are observed if γ and g 0 depart from 1, see the main panel of Fig. 18.In the latter case, we apply a phenomenological description of the density of states using the generalized normal distribution, where E is the energy and Γ is the Gamma function, while µ and σ are the mean and the standard deviation of the energy distribution, respectively.We determine the parameter β numerically and it controls the peakedness of the distribution with respect to its tails.For β = 2, one obtains the standard normal distribution, while the limiting case for β → ∞ is the uniform distribution.The variance σ 2 of the energy distribution can be estimated analytically as suggesting that at α < 1, it is extensive (Γ 2 ∝ L) by the virtue of the last term in Eq. ( 1).The density of states ρ(E) in the UM is for different parameter regimes shown in Fig. 19.Results at N = 1 in Fig. 19(a) show that ρ(E) is well described by the generalized normal distribution ρ gen (E) from Eq. (A1) with β ≳ 3, i.e., it is not accurately described by the Gaussian distribution, which would correspond to β = 2.The typ [panels (a-c)] and their derivatives dS C2) which suggest that S (q→0) → 1, as long as λ 1 is sufficiently away from 1.This is expected to be the case in the ergodic phase and in the vicinity of the critical point.In the limit q → ∞, one can simplify ⟨⟨ln(λ q 1 +(1−λ 1 ) q )⟩⟩ in Eq. (C1) by replacing it with q⟨⟨ln λ 1 ⟩⟩, which is a reasonable approximation if λ 1 is sufficiently larger than 1/2.This is a natural assumption for the nonergodic phase and, as shown in Sec.III D, also in the vicinity of the critical point.It then follows that the leading term at q → ∞ is S (q→∞) ≈ − log 2 λ 1 . (C3) This result suggests that, as expected, the two-level system is maximally entangled at the critical point if λ 1 ≈ 1/2, while the entanglement is vanishingly small if λ 1 → 1.In Sec.III D we showed that, at least for the model parameters under investigation, λ 1 at the critical point is neither close to 1/2 nor to 1, which gives rise to substantial, but not maximal entanglement of S (q→∞) .Results for the entanglement entropies at p = 1, shown in Figs. 8 and 9, are consistent with these considerations.
The above analysis is in Figs.21 and 22 extended to the case p = 4, i.e., to subsystems that consist of four most distant spin-1/2 particles.For both models, the results in Figs.21 and 22 suggest that the critical point can be rather accurately determined from the scale invariant point of S (q) emerging at α = α c as well as from the peak in its derivative dS (q) /dα.
Being more quantitative, we note that in the case of the UM, we actually observe no major differences between the p = 1 and p = 4 cases, since the results in Figs. 8 and 21 are virtually almost indistinguishable.In the QSM, however, some differences can be observed between the p = 1 case in Fig. 9 and the p = 4 case in Fig. 22.In particular, at p = 4 the crossing point of S (q) versus α at small system sizes N + L ≈ 10 emerges slightly away from α = α c , and it drifts towards α c upon increasing L. This effect is especially apparent at small values of q shown in Fig. 22(a).These results establish the view that the sharpest signatures of the critical point are in the QSM contained in properties of the particles that are most distant from the ergodic quantum dot.19), while the horizontal and vertical lines in (a) and (b) denote the averages S (2) and ∆ from Eqs. ( 15) and (17), respectively, at α = 0.8.In panels (c) and (d) we then extend the analysis of S defined in Eqs.(15) and (17), respectively, as averages over Hamiltonian eigenstates and different Hamiltonian realizations.
To better understand the origin of the emergence of S (2) being a well-defined function of ∆, we here study both quantities before the averages are carried out.Specifically, we define the entanglement entropy S Figure 23(a) shows both S (2) n,H versus ∆ n,H and S (2) H versus ∆ H at p = 1 and α = 0.8.As expected from Eq. ( 19), the entanglement entropy at p = 1 is a unique function of Schmidt gap already on a level of a single eigenstate.On the other hand, this is not the case at p = 4, see Fig. n,H versus ∆ n,H give rise to a wide cloud of points without any well defined functional dependence.
The next question that we then ask is at which level of averaging the second Rényi entanglement entropy becomes a well-defined function of the Schmidt gap even at p = 4, as suggested by the insets of Figs.H versus ∆ H at p = 1 and p = 4, respectively.The results suggest that at p = 4, the averaging over Hamiltonian eigenstates within a single Hamiltonian realization represents the key contribution to establishing a well-defined functional dependence of S (2) versus ∆.Still, a careful inspection of Fig. 23(d) reveals that the collapse of the data to a single function is not perfect, and hence the functional dependence of S (2) versus ∆ at p > 1 should not be considered as an exact property.

Appendix E: Further results for the fractal dimension
In Fig. 13 of the main text, we showed the results for the participation entropy S (q) typ versus the logarithm of the Fock space dimension ln D for the QSM.They were accurately described by the functional ansatz from Eq. ( 20) at the critical point as well as deep in the ergodic and nonergodic phases.Using this ansatz, we obtained the fractal dimension d (q) typ studied in Sec.IV.Analogous results for the UM are shown in Fig. 24.They are also accurately described by the functional ansatz from Eq. (20), i.e., the system is said to be very close to the asymptotic regime in which the L dependence of the fractal dimension d (q) typ is likely very small.The only exception may be the case for q = 0.1, see Figs. 24(a), in which S (q) typ still increases with ln D even though at larger values of q we observe no increase.

FIG. 2 .
FIG. 2. Sketch of the avalanche physics.(a) The ergodic quantum dot with N = 3 spin-1/2 particles (within a circle), which are uncoupled from the remaining L = 3 isolated spin-1/2 particles.The quantum avalanche is enabled by a hierarchical coupling of the dot to the remaining particles: (b) strong coupling to the closest particle; (c) moderate coupling to the subsequent particle; (d) weak coupling to the most distant particle.

5 FIG. 3 .
FIG. 3. Hamiltonian matrix structure of (a) the quantum sun model from Eq. (1), and (b) the ultrametric model from Eq. (4).In both cases, we consider N = L = 3 and the coupling parameter α = 0.85, while the other parameters of the Hamiltonians are g0 = γ = 1 in (a) and J = 1 in (b).A common feature of both models is the existence of 2 L = 8 dense blocks of size 2 N × 2 N = 8 × 8, denoting the ergodic quantum dot.Outside these blocks, the matrix is sparse in (a) and dense in (b).The matrix elements in each dense block in (a) are identical, except for their diagonal values, which are related to the distributions of random fields hj, see Eq. (1).The matrix elements of each block in (b) are sampled independently.When the size of the blocks is increased, the magnitude of their matrix elements is reduced accordingly.

e r g o d i c m u lt if r a c t aFIG. 15 .
FIG. 15.(a) Fractal dimension d

Figure 14
Figure 14 actually shows d (q)

14 FIG. 16 . 3 FIG. 17 .
FIG.16.The average gap ratio r defined in Eq. (9) in the UM at N = 1, as a function of α and for different system sizes L. (a) J = 0.1, (b) J = 0.5, (c) J = 1.The vertical dashed lines denote the prediction for the critical point α = αc from Eq. (6).The GOE and Poisson limits for r are denoted by the upper and lower horizontal dashed-dotted lines, respectively.

ACKNOWLEDGMENTS
We acknowledge discussions with I. M. Khaymovich.This work is supported by the Slovenian Research and Innovation Agency (ARIS), Research core funding numbers P1-0044, N1-0273, J1-50005, and N1-0369 (J.Š., M.H. and L.V.).W.D.R. was supported in part by the FWO through grant G098919N, by an internal KULeuven grant iBOF DOA/20/011 and by the Excellence of Science (EOS) programme (FWO and F.R.S.-FNRS) through grant EOS G0H1122N EOS 40007526 CHEQS.We gratefully acknowledge the High Performance Computing Research Infrastructure Eastern Region (HCP RIVR) consortium [120] and European High Performance Computing Joint Undertaking (EuroHPC JU)[121]  for funding this research by providing computing resources of the HPC system Vega at the Institute of Information sciences[122].

85 γFIG. 18 .
FIG. 18. Density of states ρ(E) the QSM.The inset shows results at g0 = γ = 1 and different values of α. Results are nearly indistinguishable from a normal distribution.The main panel shows the results at α = αc = 1/ √ 2, when g0 and γ are considerably away from 1.The dashed-dotted lines of the matching color show fits to the results using a generalized normal distribution ρgen(E) from Eq. (A1).The values of the β parameter of ρgen(E) are given in the legend, where values close to β = 2 indicate proximity to the normal distribution.The results shown above were averaged over N samples = 40 disorder realizations.

FIG. 23 .
FIG. 23.The role of averaging on the dependence of the Rényi entanglement entropies S (2) on the Schmidt gap ∆ in the UM at J = 1 and N = 1.(a) p = 1 and (b) p = 4, both at α = 0.8 and L = 14.The small green symbols denote S (2) n,H versus ∆n,H , see Eq. (D1), while the larger black symbols denote S (2) H versus ∆H , see Eq. (D2).The dotted line in (a) is the result from Eq. (19), while the horizontal and vertical lines in (a) and (b) denote the averages S(2) and ∆ from Eqs. (15) and(17), respectively, at α = 0.8.In panels (c) and (d) we then extend the analysis of S

( 2 )
H versus ∆H to α in the interval α ∈ [0.1, 0.98] and to different values of L, as indicated in the legend in (c).Results for most values of α are shown as gray symbols.The colored symbols specifically show the results for different system sizes at three different values of α = 0.9, α = αc and α = 0.4.The horizontal and vertical lines denote the averages S(2) and ∆ at a given L at α = αc.The results in panels (a) and (b) were obtained using averaging over N samples = 300 Hamiltonian realizations.Results in panels (c) and (d) were obtained by averaging over N samples = 1000, 400, 300 for L ≤ 12, L = 13, L = 14, respectively.
(q)n,H and the Schmidt gap ∆ n,H of a single eigenstate |n⟩ of a single Hamiltonian realization asS , ∆ n,H = λ 1 − λ 2 .(D1)Analogously, we define the entanglement entropy S (q) H and the Schmidt gap ∆ H of an average over Hamiltonian eigenstates of a single Hamiltonian realization as S In Fig.23we study the properties of the quantities defined in Eqs.(D1) and (D2) in the UM at J = 1 and N = 1.