Non-ergodic delocalized states for efficient population transfer within a narrow band of the energy landscape
Vadim N. Smelyanskiy Konstyantyn Kechedzhi Sergio Boixo Sergei V. Isakov Hartmut Neven Boris Altshuler







Google, Venice, CA 90291, USA

Google, Venice, CA 90291, USA
QuAIL, NASA Ames Research Center, Moffett Field, California 94035, USA
University Space Research Association, 615 National Ave, Mountain View, CA 94043 

Google, Venice, CA 90291, USA

Google, 8002 Zurich, Switzerland

Google, Venice, CA 90291, USA

Physics Department, Columbia University, 538 West 120th Street, New York, New York 10027, USA







We address the long-standing problem of the structure of the low-energy eigenstates and long-time coherent dynamics in quantum spin glass models.  This problem remains challenging due to the complex nature of the distribution of the tunneling matrix elements between the local minima of the energy landscape. We study the transverse field induced quantum dynamics of the following spin model: zero energy of all spin configurations except for a small fraction of spin configuration ("marked states") that form a narrow band at large negative energy.  The low energy dynamics can be described by the effective down-folded Hamiltonian that acts in the Hilbert subspace involving only the marked states. We obtain in an explicit form the heavy-tailed probability distribution of the off-diagonal matrix elements of the down-folded Hamiltonian.  This Hamiltonian is dense and belongs to the class of preferred basis Levy matrices (PBLM). Analytically solving nonlinear cavity equations for the ensemble of down-folded Hamiltonians allowed us to describe the statistical properties of the eigenstates.  In a broad interval of transverse fields, they are non-ergodic, albeit extended. It means that the band of marked states splits into a set of narrow minibands. Accordingly, the quantum evolution that starts from a particular marked state leads to a linear combination of the states belonging to a particular miniband. Analytical description of this qualitatively new type of quantum dynamics is a key result of our paper.
Based on our analysis we propose the  population transfer (PT)   algorithm:  the quantum evolution  under constant transverse field B_⊥  starts at a low-energy spin configuration and ends up in a superposition of Ω  spin configurations inside a narrow energy window.  
This algorithm crucially relies on non-ergodic nature of delocalized low energy eigenstates. In the considered model the runtime of the best  classical algorithm (exhaustive search) is   t_=2^n/Ω. For   √(n)≫ B_⊥≫ 1, the   typical runtime of the quantum PT algorithm  √(t_) e^n/(2B_⊥^2) scales with n and Ω as that  of  the  Grover's quantum search, except for the small correction to the exponent.  
Unlike the Hamiltonians  proposed for analog quantum unstructured search  algorithms, the model we consider
is non-integrable and the transverse field delocalizes the marked states.
As a result, our PT protocol does not require fine-tuning of the transverse field  and may be initialized in  a
computational basis state.  We find that the runtimes of the PT algorithm  are distributed according to  the  alpha-stable Levy law  with tail index 1. We argue that our approach can be applied to study PT protocol in other transverse field spin glass models, with the potential quantum advantage over classical algorithms.

Non-ergodic delocalized states for efficient population transfer within a narrow band of the energy landscape
    Boris Altshuler
    December 30, 2023
=============================================================================================================








§ INTRODUCTION 



The idea to use quantum computers for the solution of search and discreet optimization problems has been actively pursued for decades, mostly notably in connection to Grover's algorithm <cit.>,  quantum annealing
<cit.>,
and more  recently, quantum approximate optimization
<cit.>.  Quantum tunneling of collective spin
excitations was  proposed and studied experimentally as  a mechanism
for moving between states in the energy landscape  that can lead to
shorter transition time scales  compared to classical Simulated
Annealing approaches under certain conditions <cit.>. Experimental evidence of the faster time scales was later corroborated numerically using an imaginary-time Quantum Monte Carlo (QMC) algorithm <cit.>.

 Furthermore,  recent studies <cit.> have shown that in QMC, the tunneling  corresponds to the Kramers escape through  the free-energy barrier in an extended spin system that includes spin replicas in an imaginary time direction. As a result, the incoherent quantum tunneling rate does not have a scaling advantage over such a QMC simulation. This happens because  incoherent tunneling dynamics  corresponded to sequential transitions 
connecting individual minima, where each transition is dominated by a single tunneling
path <cit.>. 
In this paper we explore the qualitatively different tunneling dynamics where a large number of tunneling paths interfere constructively, giving rise to  "minibands" of the non-ergodic many-body  states delocalized in the computational basis (i.e. in the Fock  space). We demonstrate that  the  transport within the minibnads can be used  for efficient quantum search in spin glass problems.


 
 



To describe the search task we start from the  binary optimization problem where the goal 
 is to find the minimum of a
classical energy function, (z), defined  over the set of 2^n
configurations of n bits (bit-strings) z=(z^1,z^2,…,z^n)
where z^k={0,1}. In quantum algorithms (z) is typically encoded in an
n-qubit Hamiltonian

    H_ = ∑_z (z) |z⟩⟨z|

diagonal in the basis of states |z⟩  called the computational basis. Hard optimization problems have  their counterparts  in spin glass models of statistical physics <cit.>.
The energy  function  of a  hard optimization problem is characterized by a large number of spurious local minima. Low-energy minima can be  separated by a large Hamming distance (number of bit flips  transforming  one to another). Such landscape gives rise to an interesting computational primitive: given an initial bit-string z_j with atypically low energy, we wish to produce other bit-strings with energies in a narrow range Δ_ around the initial one. In general, this can be a difficult search problem if the number of bit-strings of interest is exponentially small compared to 2^n. 

Inspired by the Hamiltonian-based approaches to quantum search <cit.> and optimization  <cit.>
 we propose the following quantum population transfer (PT) protocol: first preparing the system in a computational state |z_j⟩ with classical energy (z_j), we then evolve it under the Hamiltonian

    H = H_ +H_D,   H_D =- ∑_k=0^n σ_x^k ,

without fine-tuning the evolution time nor the strength of the
time-independent transverse field .  At the final moment we
projectively measure in the computational basis and check if the
outcome z is a “solution”, i.e.,  z  z_j and the energy
(z) is inside the window Δ_. The second term in
the Hamiltonian  (<ref>) proportional to  is responsible for
the PT. It is usually referred to as a “driver Hamiltonian” in
the Quantum Annealing literature <cit.>.


We note that the  output of PT  z can be used as an input  of a classical optimization heuristic such as simulated annealing  or parallel tempering in a “hybrid” optimization algorithm <cit.> where quantum and classical steps  can be used sequentially to gain the complementary
advantages of both  <cit.>.

For random optimization problems  diagonal matrix elements (z) of the Hamiltonian (<ref>) correspond to  a problem  instance  sampled from a particular statistical ensemble. 
Since off-diagonal matrix elements  
connect states separated by one bit-flip,  Eq. (<ref>) describes the Hamiltonian of the tight-binding model
with  diagonal disorder. The underlying lattice for this model is Boolean hypercube <cit.> where individual sites correspond to bit-strings. The  model (<ref>) can be viewed as  a generalization of the Anderson model  initially introduced in the context of transport in finite
 dimensional lattices  <cit.>.   In this model, as well as  in the  original Anderson model,  there exist  bands of localized and extended states separated in energy by a so-called 
“mobility edge”. Originally,  extensions of Anderson model appeared in a variety of many-body problems 
 in condensed matter physics <cit.>  giving rise to the concept of many-body localization (MBL). It  was demonstrated in  Ref. <cit.>   that MBL is responsible for the failure of  Quantum Annealing to find a solution of  the  constraint satisfaction problem (although,  the detailed analysis of this effect is still  needed  <cit.>).
 
 
 
 
 
   In models of quantum spin glasses the existence of the two types of eigenstates and the mobility edge were studied in Refs. <cit.>. We expect the Andreson models on Boolean hypercube have  an intermediate phase characterized by the onset of non-ergodic delocalized states forming narrow minibands.  Such a phase has been observed in tight binding models on Random Regular <cit.> and fully connected graphs <cit.>. 



For spin glass models (<ref>) with  below the quantum spin glass transition, the probability density function (PDF) of the eigenvalues E_β of H is localized around the mean classical energy, with an exponentially decaying tail reaching towards the low energy states. 
We choose the interval of energies Δ_cl to be at the tail of the distribution, of E_β but sufficiently far from the ground state so that the typical  spacing of eigenvalues  is exponentially small  in n. Under these conditions   classical  states inside the energy window Δ_cl  are located near   deep local minima of the classical energy landscape  (z). Hamming distances  between the minima scale with n and the number of them is exponentially small compared to 2^n yet still exponentially large in n.


In this paper we apply the   PT protocol with Hamiltonian (<ref>) to a simple yet nontrivial model  of (z)  with the properties  mentioned above

    H_ = ∑_j=1^M(z_j)|z_j⟩⟨z_j| .

Here M≫ 1 marked states |z_j⟩ (n-bit-strings z_j) are chosen uniformly at random from all bit-strings of length n, with energies (z_j) 
independently distributed around -n within a narrow band of width W ≪. All other states z have energies (z)=0 and are separated by a large gap ≃ n from the very narrow band of marked states (see Fig. <ref>). This model is inspired by the impurity band model in doped semiconductors <cit.>. It also corresponds to a classical unstructured search problem with multiple marked states.



We provide a detailed  description of the  PT dynamics in the above model by developing a microscopic analytical theory of the "minibands" of non-ergodic delocalized states <cit.>. We derived  an effective downfolded Hamiltonian in the energy strip associated with the PT. Its matrix elements correspond to the tunneling between the deep local minima and obey the heavy-tailed statistics. The ensemble of downfolded Hamiltonians for PT corresponds to the preferred basis Levi matrices (PBLM).  We use the  cavity method for Levi matrices <cit.>
 to find analytically  the fractal dimension of the delocalized minibands, and  the probability distribution  of their spectral width. This allowed us to find the probability distribution and the scaling with n of the PT times.  
 
  It is crucial that the dynamics within the IB of the model (3) in the transverse field can be  non-ergodic yet delocalized in computational basis. The model is by no means unique from this point of view. We  believe that the extended but non-ergodic quantum states exist for  quantum extensions of any problem Hamiltonian, which is characterized by a classical spin glass behavior: for Random Energy Model  <cit.>, Sherrington-Kirkpatrick model <cit.>,  p-spin model <cit.>, K-Satisfiability <cit.>,  etc.  


 Indeed, the main difference between classical and quantum spin-glass models is the existence of the many-body localized (MBL) phase in the latter case. However we see no reason to expect a direct transition between the MBL and ergodic phases without intermediate non-ergodic phase similar to the case of ordinary Anderson localization in finite-dimensional space. This difference is due to the fact that the number of relevant bit-strings at a given Hamming  distance d from a given one increases for spin-glass models exponentially with d, or even quicker, whereas for finite-dimensional models this increase is only polynomial. 
 
 
 
A key challenge in developing  a theory of non-ergodic delocalized phase for  quantum spin glass models is the calculation of the statistics of the tunneling matrix elements between deep local minima separated by large Hamming distances d.  We derived  analytically its dependence on the transverse field B_⊥ and Hamming distance d using WKB theory of collective spin tunneling in asymptotic limit of large n. We demonstrated that in the delocalized phase  it is qualitatively different from that given by the leading order perturbation theory in B_⊥, known as a forward scattering approximation (FSA) that has been previously used in these  problems <cit.>. As a consequence, our results for the scaling of the PT time with n and the structure of the delocalized eigenstates in IB model are qualitatively different from the FSA predictions. 



 
In model (<ref>), the most efficient classical algorithm is purely random search with running time   ∼ 2^n. We find that  the  typical runtime of the PT  algorithm   displays the following scaling dependence  
 on  n 

    ∝2^n/ΩlogΩ^1/2e^n/(2 B_⊥^2) .

Here Ω≫ 1 is the number of  computational basis states  within the target window of energies that contribute with comparable probabilities to the quantum state  at the end of PT.
The expression applies  in the range of transverse fields n^1/2 ≫ -1=(1) (for arbitrary  see Eq. (<ref>)). 

The dependence  of  on Ω is the same as in the  multi-target  Grover quantum algorithm 
that searches for   Ω  marked states  starting from the  fully-symmetric   state  |S⟩=2^-n/2∑_z|z⟩. In the  Hamiltonian version of this algorithm <cit.>, one uses the projector to |S⟩  as a driver, H_D=w|S⟩⟨S|. This algorithm is proven to be optimal for problems without structure. We emphasize that according to Eq. (<ref>) the exponential scaling of  with n differs from that in the  Grover  algorithm by a term ∼ B_⊥^-2  that can be made arbitrary small at sufficiently large transverse fields. 


 PT algorithm is qualitatively different from  the quantum annealing, adiabatic optimization  and Hamiltonian implementation of  Grover search because it  exploits   the structure of the excited energy spectrum. The  PT Hamiltonian H (<ref>) is non-integrable and its eigenstates are delocalized in the low-energy manifold.

In  analytically tractable example considered here   the PT algorithm has  new  and potentially advantageous features compared to  the Grover algorithm whose 
Hamiltonian is integrable and all of its eigenstates but one are localized.  Therefore the  quantum  evolution resulting from the Grover Hamiltonian cannot form a massive superposition of Ω≫ 1 solutions if it starts from a computational basis state. The algorithm must always start from the state |S⟩.  Moreover, Grover's algorithm performance is exponentially sensitive to fine-tuning of the weight of the driver w on the scale δ w∼ 2^-n/2√(Ω). In contrast, the scaling of the runtime of PT  (<ref>) with n is robust to the choice of  that can take on a broad range of values for ≫1.

The nearly optimal (Grover-like) performance of the PT protocol is the consequence of the asymptomatic orthogonality between the eigenstates in the marked state subspace to the rest of the Hilbert space. This  suppresses the population transport from the marked states to the 𝒪(2^n) of  states |z⟩ with classical energies E_z=0  even at large B_⊥. Such "orthogonality catastrophe" cannot be obtained within the perturbative in B_⊥ approach such as FSA.





The paper is organized as follows. Sec. <ref> contains a qualitative discussion of the main results. In Sec. <ref> we develop a down-folding procedure to reduce the original problem to the nonlinear eigenproblem in the marked state subspace. In Sec. <ref> we calculate the off-diagonal (tunneling) matrix elements of the down-folded Hamiltonian and studied their dependence on n and Hamming distance   using Wentzel-Kramers-Brillouin (WKB) theory.
 In Sec. <ref> we develop an expansion of the nonlinear eigenproblem near the center of the IB shifted by transverse field and obtain the effective Hamiltonian   of the PT problem.

In Sec. <ref> we study the statistical ensemble of Hamiltonians . Sec. <ref> discusses numerical results. In Sec. <ref> we study the PT within the Born approximation. In Sec. <ref>
we estimate the number of states in the miniband.  In Sec. <ref> we provide an overview of the cavity method for dense random matrices.  In Sec. <ref> we solve the cavity equations and obtain the distributions of the real and imaginary parts of self-energy. In Sec. <ref> we discuss the complexity of PT  algorithm.
In Sec. <ref> we provide a comparison between PT and Grover's algorithm with multiple target states and systematic errors in oracle phase and driver weight.  In Sec. <ref> we provide a summary and concluding remarks.







§ QULATITATIVE DISCUSSION OF RESULTS




Each marked state |z_j⟩  is a deep local minimum of (z) separated from other minima by a typical Hamming distance  n/2  while the separation from the nearest market state is also extensive =O(n) for M=2^μ n and μ<1.


 The transverse field  gives rise to  multiqubit tunneling between  the states.   The tunneling amplitudes  from a given minimum to its neighbors located at a Hamming distance d decrease exponentially with d while the number of neighbors increases exponentially with d for d= O(n).  As a result, an eigenstate
|ψ_β⟩ of H associated  with the impurity band can become delocalized over a large subset  of marked states _β  with size 1≪ |_β|∝ M^α and  0<α≤ 1.  
For α=0 the eigenstate |ψ_β⟩ is localized, for α=1 the eigenstate is delocalized in the  entier space of  marked states. For 0<α<1 the eigenstate can be considered "non-ergodic" and its support  set _β is sparse in the space of the  marked states.  
We express the transition probability  from |z_j⟩ to |z⟩,

    P(t,z|z_j) =  |∑_β⟨z|ψ_β⟩⟨ψ_β|z_j⟩ e^-i E_β t |^2 ,

in terms of the eigenstates  and corresponding eigenvalues  of  H, where H|ψ_β⟩ = E_β|ψ_β⟩. In the delocalized phase, for a given state  |z_j⟩ there exists a large set 
of eigenstates |ψ_β⟩  that have peaks at  |z_j⟩.    These eigenstates  possess  important properties <cit.>: they have largely overlapping supports _β_β ≈ (z_j), and they are  close in energy thus forming a narrow mini-band.  The mini-band width Γ may be interpreted as the inverse scrambling time and determines the width of the  plateau in the  
Fourier-transform of the typical transition probability P̃(ω,z|z_j)  <cit.>.[The same plateau width  characterizes the frequency dependence of the  eigenfunction  overlap correlation coefficient K(ω)=M∑_j=1^M∑_β,β'   |⟨j|ψ_β⟩|^2   |⟨j|ψ_β'⟩|^2 δ(ω-E_β+E_β') <cit.>.]
 
 In other words, the  significant PT of P(t,z|z_j)   from the initial marked state |z_j⟩∈ into the other states of the same miniband  occurs over the time ∼ 1/Γ. The window Δ_ cl is related to the miniband width Γ. 
 

Understanding  the properties  of  non-ergodic delocalized  states  is crucial for describing the dynamics of quantum spin glasses driven by many-body coherent tunneling processes.
Developing its  microscopic theory is a challenging problem. This paper studies the  transport problem in an  "impurity band" (IB) model (<ref>) by making use of the down-folded   Hamiltonian   in  the marked state subspace derived in Secs. <ref>, <ref>.
While the original Hamiltonian (<ref>) is sparse in the basis of states |z⟩ (it couples only states  separated by  Hamming distance 1),
the   down-folded Hamiltonian  (<ref>) is a dense M× M matrix. 


The transverse field leads to a uniform shift ∼ of the marked state energies as shown in Sec. <ref>,  (<ref>)-(<ref>).
Diagonal elements of _ii are given by the marked state energies counted off from the center of the shifted impurity band. 
Their PDF  is assumed to be exponentially bounded with some width W. 


Each pair of  marked states is coupled via multi-qubit tunneling. 
The off-diagonal matrix elements  _ij=V(d_ij)cosϕ(d_ij) are completely determined by the Hamming distance d_ij between the marked states z_i and 
z_j. The amplitude V(d)  decays steeply with d, inversely proportional to a square root of nd (see Eq. (<ref>)). The phase  ϕ shown in Fig. <ref> monotonically increases  by (1) when d is changed by 1. In the analysis of spectral properties of _ij the quantity cosϕ(d_ij)  can be replaced by a random sign. The explicit form of V(d) and phi(d) is obtained using WKB theory of collective spin tunneling. At B_⊥ > 1 the tunneling paths correspond to long spin-flip sequences connecting the initial and final states. They  include many loops passing through the the states with Ez=0 that are neglected in FSA. 



 The typical matrix element between the two marked states is V_ typ∼ n^2 2^-n/2e^-n/(4). The typical matrix element between a given marked state and its nearest neighbor is also exponentially small in n but it is exponentially larger than the value V_ typ. 
This fact corresponds to  a  strong hierarchy of the off-diagonal matrix elements of _ij 
which is  a signature   of  their  heavy-tailed probability density function  <cit.>.  Such matrices are called Levi  matrices.

The  PDF of the rescaled squared amplitudes w_ij=V^2(d_ij)/V_ typ^2
  derived in Sec. <ref> is

    PDF(w)=1/w^2√(πlog w),   w∈[1,∞).

The particular form of scaling is the direct consequence of the fact that our problem has no "structure":   the tunneling matrix elements depend only on Hamming distance and marked states are chosen at random. 


The key difference  of the ensemble of matrices _ij from  Levy matrices studied in the literature <cit.> is that the dispersion,  W, of the diagonal matrix elements is much larger than the typical magnitude of the off-diagonal elements  V_ typ. Therefore _ij can be called
preferred basis Levi matrices (PBLM). 

We note that the existence of heavy tails in the PDF of the off-diagonal matrix elements of the down-folded  Hamiltonian  is due to  the infinite dimension of the Hilbert space of the original problem (<ref>) for  n→∞. This happens because the exponential decay of the matrix elements with the Hamming distance d is compensated by the exponential growth of the number of  states at the  distance d  from a given state. We believe that the  PBLM structure is a generic feature of the effective  Hamiltonians for PT at the tail of  the density of states in quantum spin glass problems. 


 
Unlike the standard Levi ensemble, the eigenstates of PBLM allow for the existence of non-ergodic delocalized states when the width W is much bigger than the largest off-diagonal matrix element in a typical row of _ij and much smaller than the  the largest off-diagonal element in a matrix 

    V_ typ M^1/2≪ W≪ V_ typ M .

For smaller dispersion W V_ typ M^1/2 the matrix eigenstates  are ergodic while for W   V_ typ M the eigenstates are localized.
Such phase diagram resembles the one in the Rosenzweig-Porter (RP)  model   <cit.>. The difference of RP from PBLM is that  the statistics of the off-diagonal matrix elements in the RP ensemble are Gaussian <cit.> rather than polynomial (<ref>). In this paper we will focus on exploring PT transfer within the non-ergodic delocalized phase, which is more likely to generalize to other models. We note that the localized phase does not support population transfer.


Because of the PBLM structure of the Hamiltonian  one can expect that the runtime of the PT protocol  will have  a heavy-tailed PDF  whose form is of practical interest. It is closely related to the PDF of the miniband widths  Γ∼ 1/.
We obtained the PDF(Γ)  by making use
of the cavity method for random symmetric matrices <cit.>. 


In previous work the
cavity equations were solved only in their linearized form, i.e.,  near the localization transition. We were able to solve  fully nonlinear cavity equations in the delocalized non-ergodic phase. We obtained the boundaries of the phase in terms of the ratio of W/V_ typ and also the form of  (Γ) inside the phase. It is given by the alpha-stable Levi distribution <cit.> with the tail index 1, most probable  value Γ_ typ= V_ typ(πΩlogΩ/4)^1/2, and characteristic dispersion 
πΓ_ typ/(4logΩ) where Ω is the typical number of states in the  miniband. This number Ω=(π M V_ typ/W)^2 is a square function of the ratio of the typical tunneling matrix element V_ typ to the level separation W/M.
In a non-ergodic delocalized phase M≫Ω≫1 and  the typical PT time ∼1/Γ_ typ obeys the condition

    (Mlog M)^-1/2≪ V_ typ∼ (ΩlogΩ)^-1/2≪ 1 .



We build on the observations made in the IB model and provide qualitative arguments that PT will have a  quadratic speed up over QMC in some quantum search problems where tunneling is a computational bottleneck.











 



§   DOWNFOLDING INTO THE SUBSPACE OF THE MARKED STATES AND  NONLINEAR EIGENPROBLEM 





The 
driver Hamiltonian H_D in Eq (<ref>)



connects bit-strings that are separated by a Hamming distance d=1.  
We note that, on one hand,  marked states  are separated by large Hamming distances d_ij with typical value d=n/2.
Therefore  a pair of marked states |i⟩ and |j⟩
is  coupled  by  elementary spin-flip processes corresponding to high
orders (H_D)^k  of the driver Hamiltonian with k ≥ d_ij. 
On the other hand, the
resolvent of the driver Hamiltonian 

    G(E)=1/E-H_D .

connects directly every pair of marked states. Furthermore,  because H_D is invariant under permutations of bits,  the  matrix elements  G_ij(E)=⟨z_i|G(E)|z_j⟩ depend only on the Hamming distance d_ij between the corresponding states. They are   exponentially small in n for extensive d_ij=(n). Therefore, one might expect that under certain
 conditions the quantum evolution stays approximately  confined to the marked state subspace and can be naturally described by the  downfolded Hamiltonian whose M× M matrix representation is dense in  the basis of marked states.

We use the identity

    G(E)H_|ψ⟩=|ψ⟩ ,

where E and |ψ⟩ is an   eigenvalue and the corresponding
eigenvector  of H. We introduce a new vector

    |⟩=√(H_)|ψ⟩

that  has no support in the subspace orthogonal to that of marked states.  Then, multiplying both parts of equation (<ref>) 
by √(H_), we obtain after simple transformations

    H_+Λ|⟩   =E|⟩,

where

    Λ =√(H_)H_D    G(E)√(H_) .

The
operator Λ  plays the role of a “driver Hamiltonian” in the
downfolded picture, and it couples states in the marked subspace.


Equation (<ref>) can be written in matrix form (see Appendix, Sec. <ref> for details)

    ∑_j=1^M_ij(E)_j=E _i,

where _i=⟨|z_i⟩ and _ij is a dense symmetric M× M matrix

    _ij(E)=δ_ij(z_i) +√((z_i) (z_j)) c_ij(E) .

Here   δ_ij  is  the Kronecker delta and 

    c_ij(E)=c(E,|z_i-z_j|)=⟨z_i|H_D 1/H_D-E|z_j⟩ ,

is a coupling coefficient that depends only on a Hamming distance 
|z_i-z_j|  between the bit-strings z_i and z_j.




We note that (<ref>) has the form of a nonlinear
eigenproblem. A solution of (<ref>) for the M-dimensional vector |A⟩  with nonzero norm requires

    [(E)-I E]=0,

where I is the identity matrix.  Because the downfolded  Hamiltonian
(E) explicitly depends on the energy E,   different roots E_β of the equation (<ref>)  correspond to different Hamiltonian matrices _ij(E_β). 





This can be understood from the fact that the original 2^n× 2^n
Hamiltonian (<ref>) couples the M dimensional marked state
subspace to the rest of the Hilbert space. Therefore, the projections
of the eigenvectors |ψ_β⟩  of  H onto  the  subspace
are, in general, neither normalized nor orthogonal. The same is true for the corresponding vectors |_β⟩=√(H_)|ψ_β⟩.
The normalization condition for the projections has the form (see Appendix <ref> for details)

    ∑_j,i=1^M1/Q_ji(E_β)ψ_β( z_ j) ψ_β(z_ i) =1.

where

    1/Q_ij(E)=√((z_i)(z_j)) d/d E(_i j(E)/E) .

This condition along with Eqs. (<ref>)-(<ref>)  completely defines the eigenvector projections onto the marked state subspace and the corresponding eigenvalues.


We observe that there are exactly M  roots E_β of (<ref>)   that   originate from M classical energies of the marked states (z_j) at =0. These eigenvalues and the corresponding eigenstates  will be the sole  focus  of our study. Here we just mention briefly that the rest  of the states  originate in the limit H_ cl→ 0 from the eigenstates  of the driver Hamiltonian whose energy levels
-(n-2m) (shown in Fig. <ref>) correspond to the total spin-x projections n-2m
∈[-n,n]. The levels -(n-2m) have  degeneracy nm,
which is partially  lifted  due to the coupling to the impurity band with M states.
The splitting of the  driver energy  levels -(n-2m)   increases  as a function of  transverse field in the vicinity of  “resonances”   with the levels
 of  the impurity band where (n-2m)≈ -n for integer values of  m. At resonance, the eigenstates of the   driver with total spin-x projection n-2m  are strongly hybridized with the marked states |z_j⟩.  
As will be discussed  below, the width of the resonances remains exponentially small in n for = O(n^0).
In Fig. <ref> we plot the evolution of the energy spectrum of the Hamiltonian H as a function of transverse field for the case of two impurity states M=2.





§  COUPLING COEFFICIENTS IN THE DOWNFOLDED HAMILTONIAN
 



The coupling coefficient  c_ij(E)≡ c(E,d_ij) for i≠ j  determines the off-diagonal matrix element of the downfolded Hamiltonian (<ref>) corresponding to  the  tunneling transition that connects marked states |z_i⟩ and |z_j⟩. In the IB model, the tunneling matrix element   depends only on the Hamming distance d_ij between the states. It can be calculated in explicit form from  Eq. (<ref>). For this we use  the basis of eigenstates |x⟩ of the driver Hamiltonian H_D|x⟩=H_D^x|x⟩ in Eq. (<ref>). They correspond to bit-strings x=(x^1,…, x^n) of individual qubits polarized   in positive x^a=0 and negative x^a=1 direction of the x axis. The eigenvalues of the driver  H_D^x=-(n-2h_x) depend only on the Hamming weight of the bit-strings x. Therefore one can perform explicitly the partial summation over  basis vectors |x⟩ in  
 (<ref>) under the conditions that  ∑_ax^a =k  for all bit positions   a such that z_j^a≠ z_i^a, and ∑_ax^a =l for all a where  z_j^a= z_i^a. Finally the result  (<ref>) can be written as a double sum over k∈(0,n-d_ij) and l∈(0,d_ij) 
 
    c_ij(E)=∑_k=0^n-d_ij∑_l=0^d_ijnkn-d_ijl(-1)^l  2^-n/1+E/(n-2 k-2 l).

 Here d_ij is the Hamming distance between bit-strings  z_i and z_j. 
 Plots of coupling coefficients as a function of Hamming distance d based on (<ref>) are given in Fig. <ref>. They display qualitatively different behavior  depending on the value of the parameter n/|E|.
 
 
 
 For n/|E|<1 the coefficient c(E,d) decays exponentially  with d in the entire range of values d∈[0,n].  For  n/|E|>1 the coefficient decays until d∼ n/2, corresponding to minimum overlap between the marked states, and then begins to grow. For large transverse field ≫1 the behavior  with d is nearly symmetric  with respect to d=n/2 and to leading order it does not depend  on .  Unfortunately, the  expression  (<ref>) is  quite involved and is not suitable for the  study of the  asymptotic properties of the  population transfer in the  limit of large n.

For a very weak transverse field ≪ n^-1/2 using perturbation theory in  to the leading order  one can obtain a standard expression   <cit.> for the coupling coefficient, |c(E,d)|≃ d ! (/n)^d.   It is given by the sum of the transition amplitudes over the  d! shortest paths between the states |z_i⟩ and |z_j⟩ separated by a Hamming distance d. Intermediate states |z⟩ along each path correspond to (z)=0 while energies of initial and final states are -n (see Fig. <ref>).


For larger transverse field values  (but still ≪1) the perturbative expression in the small- limit can be modified to include the range of =(n^0) but ≪1.
In that range 

    |c(E,d)|≃ d !/n^d e^B_⊥^2d^2(3n-2d)/6n^2+n/12 .

One can see that  for small   matrix element falls down with d extremely steeply despite the presence of the 
factorial factor d ! in (<ref>). We note that this  perturbation (FSA) expression is qualitatively valid  in the range B_⊥<|E/n|≪ 1. It gives a correct leading order form of the mobility edge in quantum REM <cit.> at small B_⊥≪ 1. 

For transverse field, >|E|/n, the dependence  of c_ij(E) on d_ij changes qualitatively.  It becomes non-monotonic, reaching its  minimum at the point n/2 of minimum overlap between the bit-strings z_i and  z_j.  In a certain region around the minimum  it has oscillatory behavior, as seen in Fig. <ref>. The boundary of this region is shown with black dots. The details of the behavior in the oscillatory region are shown in  Fig. <ref>. 
The exponential dependence of  the envelope of c(E,d) on d is captured by  the factor 1/nd and is  independent on the transverse field strength.  This region of d and values of >|E|/n are of the most relevance to the transport in non-ergodic minibands which is of central interest to in this paper.



 §.§ WKB calculation of coupling coefficients


 In this paper we   develop an approach (described in the Appendix <ref>)  based on the   WKB theory for large spin  <cit.> to calculate the coefficient c(E,d) for n≫1 and arbitrary values of transverse  fields   without relying on perturbation theory in .
The coefficient  c_i j(E) can be expressed in terms of the operator of the total  spin-x projection S_x=1/2∑_j=1^nσ_x^j

    c_i j(E)=δ_ij - E ⟨z_i|( E+2 S_x)^-1|z_j⟩ .



We will utilize   the basis of  eigenstates  |m⟩ of the operator S_z=∑_k=1^nσ_z^k corresponding to  its eigenvalues   m∈[-n/2,n/2] and  the maximum value of the total spin S=n/2

    S_z |m⟩ = m|m⟩,    m=-n/2 ,…,n/2 .

The  state  | n/2-d⟩ 
 is a  normalized sum  of all computational basis states |z⟩ with  d spins pointing in the  negative z direction and n-d spins   pointing in the positive z direction (m=n/2-d)
 
    | n/2 -d⟩=1/√(nd)∑_z∈δ_|z|, d|z⟩.

Here |z|=∑_k=1^nz^k and δ_k,d is a Kronecker delta.


Because the coefficients c_i j(E) (<ref>) depend only on the Hamming distance |z_i-z_j| between the bit-strings z_i and z_j, we can assume,  without loss of generality, that  in  Eq. (<ref>) one of the bit-strings, e.g., |z_j⟩,  corresponds to all individual spins  pointing in the positive z direction

    |z_j⟩=|00 … 0⟩≡|n/2⟩  (m_z=n/2)

The main observation is that  we can  pick, instead of the state |z_i⟩, any computational basis state |z⟩ whose Hamming weight satisfies the condition |z|=|z_i| without changing the value of the coefficient c_ij(E)=c(E,|z_i|).  Therefore averaging both  sides  of the Eq. (<ref>) over the states |z_i⟩ that satisfy the condition |z_i|=d for some integer d∈[0,n]  we obtain





    c(E,d)=δ_d,0-E/√(nd) G_n/2 -d,n/2(E)     .

Here G_m,n/2(E)=⟨m|(E+2 S_x)^-1|n/2⟩  are the  matrix elements of the resolvent  (<ref>) of the transverse field Hamiltonian  H_D between the states (<ref>) that belong to  a maximum total spin subspace S=n/2.



 


As will be shown below, for typical instances  of the ensemble of Hamiltonians H, the Hamming distance from  a randomly selected marked state to its closest neighbor is an extensive quantity (n). Therefore the  above off-diagonal  matrix elements of the resolvent can be analyzed in a semiclassical approximation  corresponding to  S=n/2 ≫1. This  approximation for the quantum propagator of a large spin and diagonal elements of the resolvent was considered in <cit.> using the spin coherent state  path-integral representation. The analysis in these papers was   quite involved because the path-integral formulation requires a careful treatment of the fluctuation determinant and a so-called   Solari-Kochetov correction in the action. Also, these results were  focused on  a general case of  large spin Hamiltonian
 and  only considered  diagonal elements of the resolvent. Because of this,  instead of trying to extend the results in <cit.>  to our case, we follow a different path.
 
The resolvent satisfies the equation 
    I-2 S_x G(E)=E G(E)
 where I is the identity operator. We write this equation in the basis of states |m⟩    (<ref>). From (<ref>)  we obtain

    δ_m,n/2+∑_s=± 1u(m-s/2)  G_m+s,n/2=E G_m,n/2,




    u(m)=-√(L^2-m^2),   L=n+1/2.

In the limit of large n≫1 we  solve this equation using the discrete Wentzel-Kramers-Brillouin (WKB) approximation  method <cit.>. In the WKB analysis of Eq. (<ref>) the function  2u(m) plays the role of an effective potential for the classical system with coordinate m and energy E. For 2u(m)>E the WKB solution for the resolvent G_m,n/2(E)  displays an oscillatory behavior with m while for 2u(m)<E it exponentially increases with m. The boundaries of the oscillatory region m∈[-m_0(E),m_0(E)] are “turning points” of the classical motion and  are given by the condition 2u(m_0)=E (see Fig. <ref>) where 

    m_0=√(L^2-E/4^2) .

In Fig. <ref> we plot the comparison between the coefficient c(E,d) computed based on the exact expression (<ref>) and the WKB asymptotic  (details of the WKB analysis are given in Appendix <ref>).

  


In what follows we will be interested in the region d∈[n/2-m_0,n/2+m_0] with  m_0≃√((n/2)^2-(E/)^2)  defined by the condition u(m_0)=E.
This is the region of oscillatory behavior of c(E,d) with d where  the  leading order exponential dependence on n   is given by the expression

    c(E,d) ∝1/√(nd)  e^-nθ()sinϕ(E,d) ,

with the prefactor  given in Appendix, Eqs. (<ref>),(<ref>).


The function θ() in (<ref>) equals

    θ() =2arctanh(B_⊥^-1)+ln1-B_⊥^-2/4 .

It behaves  at large argument  as θ≃ 1/(4).

 




An explicit form of  the WKB phase ϕ(E,d) in (<ref>) is given in Appendix, Eq. (<ref>). The dependence of the phase on d  for different values of  is shown in the Fig. <ref>. 
This phase varies by (1)  when d is changed by 1 and it is  responsible for fast oscillation of the coupling coefficient with the Hamming distance between marked states d. Its dependence on d simplifies in the limit of large transverse field ≫1:

    ϕ(E,d)≃π d/2-π n/4χ(1/2-d/n)/ ,

where  χ(x) ≃ 1-2 arcsin(x)/π+(n^-1). The  second term  in (<ref>)  is much smaller the the first one, and  varies very little when d is changed by 1. A predominately linear dependence of ϕ(E,d) on d at large fields can be seen in Fig. <ref>.
This property will be important in the analytical study of population transfer.







For large transverse fields  the  magnitude of  the squared coupling coefficient  (<ref>)  can be estimated to exponential accuracy as c^2(E,d)∼exp[-n/(2B_⊥^2)] /nd. 
We note that the number of marked states M_d accessible via all possible d-bit flips from a given state is M_d=M 2^-nnd. Therefore  the leading order dependence of the coupling coefficient on d is 
 proportional to  1/√(M_d). As will be shown later,  in the limit of large transverse fields this leads to a nearly Grover complexity of the PT algorithm, up to a factor 
∼exp[-n/(4B_⊥^2)], which gives very small correction to Grover scaling for large B_⊥. However  when d decreases below the boundary value
d<n/2-m_0,  the coupling coefficient grows exponentially faster than 1/√(M_d), as follows from the discussion in Appendix (cf. Eq. (<ref>)). 

 







  §.§.§ Coupling coefficients in the region  | n/2-d  | ≪  n



In this region |m|≪ n and we can write

    ∫_m^m_0(E) p(k,E)dk    -π/4≃ϕ_∞(E)
        -marcsin(2m_0(E)/n)

Plugging this into (<ref>), (<ref>) we obtain in leading order in n
 
    f_m=f_0 sin(ϕ_∞(E)-marcsin(√(1-B_⊥^-2))/sinϕ_∞(E)

 We consider the limit
 
    ≫ 1   (= O(n^0)

 In this limit 
 
    sin(ϕ_∞(E)-marcsin(2m_0(E)/n))≃

 

 
    ≃ [[ (-1)^m/2sin(ϕ_∞(E)+m/),                 m- even; (-1)^m+1/2cos(ϕ_∞(E)+m/                  m- odd ].

 With this the expression for the squared coupling coefficients c^2(E,d) has the form
 
    c^2(E,d)= A(E) (1-(-1)^m cos(2ϕ_∞(E)+2m/))/nd.

 


 

§  DOWNFOLDED  HAMILTONIAN NEAR THE  CENTER OF THE IMPURITY BAND 

  
 
The coupling coefficients  c(E,d) (<ref>) decay exponentially with Hamming distances for d= O(n) (see details in Sec. <ref>).
Marked states are selected at random and Hamming distances between them are order n when the number of the states M is exponentially smaller than 2^n. 
Because the off-diagonal matrix elements of the downfolded  Hamiltonian H_ij(E)∝ c(E,d_ij) they are exponentially small in n.
At the same time  the width  of the distribution of energies  of the marked states (z_j)=-n +ϵ_j  is also assumed to be very small, W≪  (it is exponentially small in n for the cases of interest).
Therefore we  can solve the nonlinear eigenvalue problem  (<ref>)-(<ref>)  by an iterative approach  
treating the off-diagonal part of H(E)  and terms ∝ϵ_j as a perturbation. Details are given in Appendix <ref>.




At zeroth-order in the perturbation, the down-folded Hamiltonian  _ij^(0)(E)= δ_ij n(c(E,0)-1) has one M-fold degenerate
energy level E^(0) that is a root of the equation  _ij^(0)(E)=E that originates from the   marked  state energy, E^(0)→ -n,
in the limit of → 0. 

Using c(E,0) from Eqs. (<ref>), (<ref>) the explicit form of the equation  for E^(0) is

    E^(0) =-n-Δ_0 ,



 
    Δ_0=n 2^-n∑_d=0^nnd(n-2d)/n+Δ_0-(n-2d)  .
 
 Here    Δ_0 is  the root of the  above  transcendental equation that satisfies the condition  lim_→ 0Δ_0=0. In general,  
  the sum (<ref>) is dominated by the region of   values of d such that |d-n/2|= O(n^1/2) where the factor  2^-nnd  
reaches its maximum ∼ n^-1/2. In that region we replace the binomial coefficient with a Gaussian function of d and  the summation with the integral over d. Taking the integral 
we    obtain   Δ_0  in a form of a series expansion in powers of n^-1 

    Δ_0≃ -B_⊥^2-B_⊥^4/n+ O(n^-2),
  
A comparison between the exact and asymptotic solutions for  Δ_0 is shown in Fig. <ref>. For ≪ n^1/2 the overall shift of the energies of the marked states is negative and quadratic in .





According to Eq. (<ref>)   all  M  degenerate eigenstates  |ψ⟩_β have the same weight Q(E^(0))=∑_j=1^M|ψ_β^(0)(z_j)|^2 on the marked state subspace.   In the  large n limit we have

    Q(E^(0))≃ 1-B_⊥^2/n+ O(B_⊥^4/n^2).

Under the condition   ∝ B_⊥^2/n≪1, the  eigenstates  are dominated by their projections on the marked state subspace. In the limit n→∞ they are 
 asymptotically orthogonal to the computational  basis states outside the IB.  Such "orthogonality catastrophe" cannot be obtained within the perturbative in B_⊥ approach such as FSA. 
 
 The exact dependence of the weight Q on transverse field  is given in Fig. <ref>. The  expression (<ref>)
  is valid for  away from their "resonant" values B_⊥,p≃ n/(n-2p) where the M-fold degenerate energy level   "crosses”  the eigenvalues of the driver  Hamiltonian, E^(0)=-(n-2p), for integer values of p, as shown in Fig. <ref>. The width of such resonance regions Δ B_⊥,p∝ 2^-n/2np  
  remains exponentially small in n  for   n/2-p≫ n^1/2. 
  
  




 In this study we will  focus on the off-resonance case depicted in Fig. <ref>. One can  see from  Fig. <ref> that B_⊥,p increases with p and so is the width of resonance region.  For   parametrically large  compared to unity one needs to make sure that n is also large enough so that the width of the resonance regions  is small (cf.  Fig. <ref>). Away from resonance, all M impurity band eigenstates are  well localized in the marked states subspace (cf.  (<ref>)).
 
 





 §.§  Vicinity of resonances


We observe  that when  is increasing, E^(0)  repeatedly “crosses” one of the eigenvalues of the driver  Hamiltonian, E^(0)=-(n-2m) for integer values of m, as shown in Fig. <ref>. As a result the denominator in Eq. (<ref>) corresponding to d=m will become zero at a “resonant” transverse field value =B_⊥  m  which is a root of the equation

    B_⊥  m=n+Δ_0(B_⊥  m)/n-2m,

approximately given by B_⊥  m≃ n/(n-2m).  For the magnitude of the transverse fields ≪√(n) where  perturbation theory applies, levels of the driver that can cross with the impurity band satisfy the condition n/2-m≫ n^1/2.

Near the mth resonance the term with d=m in the sum (<ref>) becomes anomalously large due to a small denominator despite the factor p_m being very small. One can show that the width of the resonance region is given by

    Δ B_⊥  m∼   2^-n/2nm B_⊥  m .

There the weight factor Q is decreasing dramatically (cf.   Fig. <ref>) and   the above perturbation theory breaks down. This corresponds to   the hybridization of the marked states  with the eigenstates of the transverse field operator H_D corresponding to its  eigenvalue - (n-2m).   The width of the resonant regions  Δ B_⊥  m   (<ref>)  grows   rapidly with m due to the increase in the degeneracy of the corresponding eigenvalue of  H_D.  However  it   remains exponentially small in n  for   n/2-m≫ n^1/2.

 





In the spirit of the degenerate perturbation theory, there exists an effective Hamiltonian  that  determines the correct zeroth order eigenstates 
and removes the degeneracy of the energy levels

    |ψ_β^(0)⟩=E_β^(1)|ψ^(0)_β⟩ ,

Its matrix in the basis of the marked states has the form _ij=δ_ij_i+n c(E^(0),d_ij) where we neglected small non-important corrections (see Appendix <ref>). Using the  expression for  the coupling coefficient 
(<ref>)  given in Appendix <ref>, ((<ref>),(<ref>)) we have

    _ij=δ_ijϵ_j+(1-δ_ij)_ij√(2)sinϕ(d_ij) .
 




Here   ϕ(d)≡ϕ(E^(0),d) is a WKB phase shown in Fig. <ref> that describes the  oscillation of the matrix elements with the Hamming distance.
Its explicit form is given in Appendix <ref>, Eq. (<ref>) and also above in Eq. (<ref>) for the case of large transverse fields. The amplitude _ij
equals 

    _ij≡ V(d_ij),   V(d)=√(A(d/n))   n^5/4   e^-nθ()/√(nd) ,

 where i ≠j and the coefficient A(ρ) equals (cf. 
  (<ref>))
  
    A(ρ)=√(π/32 )e^-/(B_⊥^2-1)υ(ρ)sin^4(ϕ(n/2)) ,


    υ(ρ)=(1-(1-2ρ)^2/1-B_⊥^-2)^1/2 .

It is independent on n apart from the phase
  ϕ(n/2) whose  explicit form is

    ϕ(n/2)=π/4(n(1-B_⊥^-1)-) .

 The function θ() is given in (<ref>).
Expanding  (<ref>) in the limit  ≫1, 

    θ≃1/4 B_⊥^2+1/24 B_⊥^4+1/60B_⊥^6+….

In that limit θ≪ 1. We note that even for modest values of transverse field, e.g., ≃ 1.46 (corresponding to that in the Fig. <ref>) the first term provides a good estimate to the value of θ≃ 0.13 (error 9%). We shall refer to  in (<ref>) as the Impurity Band (IB) Hamiltonian.


The form of the IB Hamiltonian (<ref>) only applies to  the region of oscillatory behavior  d_ij∈[n/2-m_0,n/2+m_0]  of the coupling coefficients c_ij(E) with Hamming distance d_ij where  m_0  is given in (<ref>). This above condition for d_ij is always satisfied   in a typical row of the matrix d_ij for the  values of M  considered in the paper (see the discussion in Appendix <ref> and  Eq. (<ref>)).







§  STATISTICAL ENSEMBLE OF THE  IMPURITY BAND HAMILTONIANS
 

Properties of the eigenstates and eigenvalues (<ref>) of the IB Hamiltonian  (<ref>)
  determine the population transfer within the Impurity Band and are thus 
  of the central interest for us in this study.  They depend on the statistical ensemble of IB Hamiltonians. In the model  considered in this paper diagonal elements  ϵ_j of    are  selected at random, independently  from each other and from  the choice of the corresponding marked states  |z_j⟩.    In the present discussion we assume
  that  the   PDF 
    p()  of  ϵ_j   is exponential bounded with the width  W.   The results do not depend on the particular form  of p(). For the sake of specificity in calculations we  will use  the window function form 
    
    p()=θ(W/2-||) ,

where θ(x) is a Heaviside theta function.
For the physical effects discussed in the paper to take place the width W needs to scale down exponentially with n  

    lim_n→∞log (W^1/n)= O(n^0).








 §.§ Off-diagonal matrix elements




For fixed energies ϵ_j the  matrix of the IB Hamiltonian _ij is entirely determined by the symmetric matrix of Hamming distances d_ij between the bit-strings corresponding to the marked  states. The set of M bit-strings is randomly sampled from the full set of all possible  2^n bit-strings  without replacement, see Appendix <ref>. Elements of the matrix d_ij above or below the main diagonal will be considered independent from each other and taken from  the binomial distribution p_d,

    p_d =1/Z 2^-nnd,   Z=∑_d=1^n2^-nnd,
 
under condition 1≪ M ≪ 2^n/2. Then, for a given row of the matrix M× M of Hamming distances d_ij,   the numbers of elements M_j^(d) with d_ij=d are samples from the multinomial distribution with mean values ⟨ M_j^(d)⟩=M p_d (see Fig. <ref>).   According to  (<ref>),(<ref>) the statistical ensemble of  IB Hamiltonians (<ref>) corresponds to that of symmetric random matrices whose associated graphs are fully connected and  matrix elements are  statistically independent.  




As will be seen below  the spectral properties of   that are relevant for  our study  are determined by _ij^2 and not by the oscillatory factor in (<ref>). Therefore we will be interested in the PDF of  _ij^2 

    P(_ij^2)= ∑_d=1^n p_d   δ(V^2(d)-_ij^2) ,

where i≠ j.

 


  §.§.§  Typical and extreme values of the off-diagonal matrix elements _ij

  For a  randomly chosen row of the  matrix of Hamming distances d_ij the most probable value (mean)  of its elements equals to n/2.   According to  (<ref>),  the    off-diagonal matrix elements   _ij
  decrease rapidly with the Hamming distance d_ij, reaching the  minimum  value at d_ij≃ n/2. Therefore a typical minimum value of the matrix elements _ij  corresponds to  a typical value overall. We estimate it using Eq. (<ref>)  and Stirling's approximation 
  
    V_ typ=V(n/2)≃(π A^2  /2)^1/4  n^2 2^-n/2 e^-n θ .

 where coefficient A=A(E^(0),1/2)   (<ref>) is essentially n-independent between the resonances  and θ is given in (<ref>). The matrix elements  _ij that  scale with n as the typical value in (<ref>)  correspond to  | n/2-d_ij| = O(√(n)).  
 

 We note that in the Fig. <ref> the plot points do not reach the boundaries of the interval d=0,n. In the matrix of Hamming distances   d_ij the  typical  smallest off-diagonal  element in a  randomly chosen row can be estimated  as follows M p_=1 where p_d is binomial distribution (<ref>)

    min_j≠ i, 1≤ j≤ M d_ij∼,    M 2^-nn =1.

  Using Stirling's approximation for  factorials in the limit n≫ 1 it is easy to show that minimum Hamming distance in a row is extensive for M=2^μ n, μ <1. 
 
 The typical largest magnitude off-diagonal matrix element in a randomly chosen row of _ij is equal to V().  Using Stirling's approximation in (<ref>) we get,
 
    max_j≠ i, 1≤ j≤ M|_ij| ∼  M^1/2  V_ typ .

 Using (<ref>) one can see   that the maximum off-diagonal matrix element in a randomly chosen row is still exponentially small in n.
 
 Similarly, one can estimate  the typical value of the absolute minimum  of  a Hamming distance d_ij  between a pair of marked states. This distance remains extensive for μ<1, M=2^μ n. This distance corresponds to the overall largest in magnitude element  of the matrix _ij 
 
    max_1≤ i< j≤ M |_ij|  ∼  M  V_ typ .

 Using (<ref>) the largest element is exponentially small  in n provided that   μ< 1/2 which corresponds to the condition of statistical independence of the elements of _ij. A tight bound for the maximum eigenvalues of H can be obtained using Gerschgorin circle theorem <cit.>, see Appendix <ref>.


 
 
 

 §.§ Heavy tails

 
 It can be shown that the variance of _ij  is not a good
statistical characteristic of its PDF and is  dominated by the extremely rare atypical instances of the
ensemble (see details in the  Appendix <ref>). We observe that the relationship between the typical matrix element (<ref>), maximum matrix element in a randomly chosen row of _ij (<ref>), and the  largest element of _ij overall (<ref>) form a strong hierarchy that is a characteristic of the ensemble of    dense   matrices with broad non-exponential distribution of matrix elements (Levy matrices) <cit.>. The form of  the hierarchy <cit.>  suggests  (up to a logarithmic factors)  the following asymptotic behavior at the tail  of the PDF of the matrix elements: 
 
    (_ij^2)∝ |_ij|^-2,

  for |_ij|≫ V_ typ.
 

 We will build on the above observation and  obtain the explicit form of the  PDF of the matrix elements P(_ij^2) (<ref>),  including its tails. 
 In the asymptotic limit of large n≫ 1 we consider n to be a continuous variable (the validity of this approximation will be justified below). We replace the summation over d in (<ref>) by an integral and Kronecker delta 
 δ(x) by Dirac delta
 
    P(_ij^2)≃∫_0^n p_x   δ(V^2(x)-_ij^2)dx .

 This expression is obtained  using the analytical continuation  of the binomial distribution  p_d (<ref>)  from the integer domain d∈ (0,n)  onto  the interval of a real axis x∈ (0,n) in terms of the Beta function and the resulting identity  ∫_0^ndx  p_x=1 (see Appendix <ref> for details).  


In what follows we will study the   rescaled quantities 

    w_ij=V_ij^2/V_ typ^2≡V^2(d_ij)/V_ typ^2 ,

where   i≠ j and  V_ typ  is  given in (<ref>). We apply Stirling's approximation for  the binomial coefficient in Eq. (<ref>) and (<ref>) and obtain asymptotic expressions for   V^2(d) and p_d, respectively. Plugging them into the  (<ref>) and taking the integral there we  can obtain   the PDF

    g(w_ij)=V_ typ^2 P(V_ typ^2 w_ij) .

whose form is given in Appendix,  Eqs. (<ref>),(<ref>).


The  following assumption will be applied throughout the paper

    M=2^μ n, μ≪ 1.

According to Eqs. (<ref>), (<ref>),(<ref>)
a typical largest  element in a randomly chosen row of the matrix w_ij is ∼ M. Therefore based on (<ref>) the following condition is satisfied in a randomly chosen row of w_ij

    1/nlog_2 w_ij≪ 1   (1≤ w_ij M) .

Under this condition,  the PDF of w_ij  takes  a particularly simple  form, g(w)≃ g_∞(w)  
 
    g_∞(w)=1/w^2√(πlog w),   w∈(1,∞) ,

with normalization condition ∫_1^∞g_∞(w)dw=1. Details of the derivation are given in Appendix <ref>. 

The  above analysis assumes the   scaling behavior (<ref>) of  _ij  with d_ij that requires  |n/2-d_ij| <m_0  with m_0 given in (<ref>). As shown in  Appendix <ref>  this condition is always satisfied for a typical row of d_ij provided the constraint (<ref>) on the values of M. 
 







 §.§ Preferred basis Levy matrices (PBLMs)


The problem of population transfer is reduced to the analysis of the described above ensemble of real symmetric M× M matrices  _ij of the down-folded IB Hamiltonian (<ref>). 

 The matrices _ij form an ensemble of preferred basis Levy matrices (PBLMs), a generalization of Levy matrices actively studied in the literature (cf., e.g., <cit.>). Unlike Levy matrices PBLMs have a new control parameter: the ratio of typical diagonal to off-diagonal matrix elements W/ that controls the preferential basis (computation basis). This distinction is analogous to that between Gaussian Orthogonal ensemble (GOE) and the Gaussian ensemble with broken SU(N) symmetry, the Rosenzweig-Porter (RP) model <cit.>.

Recent studies of RP ensemble <cit.>  demonstrated two localization transitions that occur with varying parameter that controls the relative weight of the diagonal and off-diagonal matrix elements. One of them is the Anderson transition from localized to the extended states that are non-ergodic and posses distinct multifractal features. These states and the corresponding eigenvalues are organized in  "minibands” so that the states within the same miniband mostly share the same support over basis states. The spectral width of the minibands is polynomially  small (in M) compared to W.  The second transition is from the extended non-ergodic states to the extended ergodic states similar to the eigenstates of the Gaussian Orthogonal Ensemble. We demonstrate analogous behavior in the IB model and analyze the population transfer in the non-ergodic regime. 




§  NUMERICAL SIMULATIONS:  MINIBANDS OF NON-ERGODIC DELOCALIZED STATES


In this Section we report exact diagonalization analysis of both the eigenvector statistics and the dynamical eigenstate correlator. Instead of the sparse 2^n×2^n Hamiltonian Eq. (<ref>), it is efficient to diagonalize the dense M×  M matrices obtained by down-folding the Hamiltonian into the marked states subspace. This allows access to systems of n=200 qubits, reducing the finite size effects. The down-folded matrix Hamiltonian ensemble, is constructed as in Sec. <ref>, 

    _ii=ϵ_i,   _ij=n c(E^(0),d_ij),

where the diagonal elements ϵ_m are distributed uniformly in the energy window [-n-W/2, -n+W/2 ], and the off-diagonal elements are constructed by sampling Hamming distances between uniformly random bitstrings of length n and using Eq. (<ref>) with E=E^(0) determined from Eqs. (<ref>),(<ref>). 

We introduce the scaling of the width of the distribution of ϵ_m with the matrix size M,

    W= λ M^γ/2 V_typ ,

where γ is a real non-negative parameter that controls the scaling of the typical diagonal to off-diagonal matrix element V_typ given in Eq. (<ref>), and λ is an auxiliary constant of order one.




 §.§ Eigenvector statistics





 

We define the inverse participation ratios (IPRs) I_q and the entropy H^z as,

    I_q    = ∑_i |⟨ψ_β | i ⟩ |^2q,
    
      H^z    = -∑_i |⟨ψ_β | i ⟩ |^2 ln |⟨ψ_β | i ⟩ |^2,

where ψ_β denotes an eigenstate with eigenvalue E_β. IPR I_2 is the second moment of the wave function probability distribution |⟨ψ_β | i ⟩ |^2 in the computational basis (bitstrings) |i⟩. The entropy H^z characterizes the support set of an eigenstate in the computational basis <cit.>, i.e. the subset of bitstrings where the probabilities |⟨ψ_β | i ⟩ |^2 are concentrated.

 Fig. <ref> shows the participation ratio I_2 as a function of the ratio of mean level spacing δϵ to the typical matrix element V_ typ, a measure of the number of states in resonance with a typical classical level ϵ_i. The regime δϵ≫ V_ typ corresponds to the localized phase, where the eigenstates have significant weight on a small number of bitstrings that are close to each other in Hamming distance. In this regime  I_2 ∼ 1 and is system size independent. In our model marked states are separated by Hamming distance d≈  n/2+𝒪(√(n)) with high probability and therefore most localized states have sharp peaks at exactly one bitstring, hence I_2≈ 1. As the ratio δϵ/V_ typ decreases I_2 becomes system size dependent. Fig. <ref> indicates  that the combination I_2 M/3∼1 becomes system size independent as level spacing becomes smaller than the typical matrix element, characteristic of the delocalized regime, where the wave function amplitude spreads over 𝒪(M) bitstrings, |⟨ψ_β | i ⟩ |^2∼1/M. Saturation value of I_2 M ∼ 3 is consistent with approach to Porter-Thomas distribution of the wave function amplitudes. Both Figs. <ref> and <ref> show a wide intermediate region between the localized and ergodic phases where non-ergodic dynamics is expected. This intermediate regime becomes apparent in Fig. <ref> where we introduce the multi-fractal dimensions  D_q  and D_1 which determine the scaling of I_q and H^z with M, respectively,

    ln I_q(M)     = -D_q(q-1) ln M + c_q, 
    
      H^z(M)    = D_1 ln M + c_1,

where c_q is a q-dependent fitting parameter. The extracted dimensions shown in Fig. <ref> as a function of the parameter γ vary continuously between D_q=1 in the ergodic phase γ≤1 and D_q=0 in the localized phase γ≥2, with 1 <γ  < 2 corresponding to non-ergodic regime for q=1,2. 


 




 §.§ Eigenstate overlap correlator for non-ergodic minibands


Population transfer dynamics in the non-ergodic regime can be characterized by the survival probability, see Section <ref>. The Fourier transform of the survival probability for a given initial marked state i is given by,

    p_i( ω)=Re∫_0^∞ dt e^iω t| ⟨ i | ψ (t) ⟩|^2 
    
    =π∑_β,β'| ⟨ i | ψ_β⟩|^2| ⟨ψ_β'  | i ⟩|^2δ (E_β-E_β'-ω).

Note that the limit ω→ 0 gives the inverse participation ratio of a given bitstring in the basis of eigenstates,

    p_i(0) =π∑_β| ⟨ i | ψ_β⟩|^4.
  
The average of p_i(ω) over the initial state is related to the overlap correlation function K(ω) defined by <cit.>,

    K(ω)    ≡1/M∑_i,β,β' |⟨ψ_β|i⟩|^2 |⟨ψ_β'|i⟩|^2δ (E_β-E_β'-ω) 
       = 1/π M∑_i p_i( ω) .

The fractal dimension extracted from the scaling of K(0) with M is shown in Fig. <ref>, it follows closely those extracted from the IPR in the computational basis. The collapse of the plots in Fig. <ref> is achieved when the frequency is rescaled by the characteristic energy,

    Γ_ε =Γ_typ M^ε,   Γ∝ V_typ M^1-γ/2 (log M)^1/2,

with a fitting parameter ε≪1. The correlator K(ω) is constant for a range of energy differences ω  <Γ_ε and decays quickly ∝ω^-2 as ω > Γ_ε. This can be interpreted in terms of the formation of non-ergodic mini-bands of eigenstates that share support in computation basis: for an average bitstring there is a range of eigenenergies E_β within a width Γ_ε around a bitstring dependent value ϵ_j where the eigenfunction overlaps with z_j are relatively large, whereas for larger energy difference the correlation decays quickly below the value corresponding to uncorrelated case K(ω)< 1/M 
 i.e. the amplitudes repel each other. 
The relation between the survival probability and eigenfunction overlap correlator, Eq. (<ref>) suggests that the characteristic population transfer is given by the inverse of the characteristic energy scale of the miniband width Γ_ε, the range of energy eigenstates with significant amplitude at the given bitstring. The auxiliary fitting parameter takes a small value ε=0.05 indicating only a small deviation from Γ_ typ most likely due to finite size effects. In Appendix <ref> we show the results of direct simulation of dynamics of the model in the course of the PT protocol and confirm the scaling of the PT time.








 §.§  Discussion of numerical results


 The size of the matrix of marked states used in exact diagonalization M≤ 20000 is a small fraction of the size of the total Hilbert space Hamiltonian 2^n×2^n with n=200. For such a small sample the distribution of Hamming distances d_ij between marked states is dominated by |d_ij-n/2|∼𝒪(√(n)). In this regime the square of the off-diagonal matrix element, see Sec. <ref>, has approximately Gaussian dependence on d_ij (cf. Eqs. (<ref>),(<ref>))

    _ij^2 ≈ V_typ^2  exp(2/n(d_ij-n/2)^2),

and the probability to find a pair of bitstrings at a smaller distance d_ij is strongly suppressed. The sign of _ij rapidly fluctuates as a function of d_ij resulting in a negligible average ⟨_ij(d) ⟩∼𝒪(2^-n). The distribution of off-diagonal matrix elements in Eq. (<ref>) is non-Gaussian and instead has a heavy tail that cannot be fully characterized by the variance alone, see Section <ref> and Appendix <ref> where we introduced the class of Preferred Basis Levy Matrices and derived the asymptotic form of the distribution of matrix elements. For numerically accessible matrix sizes M we expect the deviation from the Gaussian distribution in the observables to be very small.


The eigenstate statistics and the respective fractal dimensions for the model  Eq. (<ref>) can be calculated using strong disorder perturbation theory. The calculation proceeds similarly to that in Ref. Kravtzov2015RzPr resulting in,

    D_1=D_2=2-γ.

Comparison of the approximate Eq. (<ref>) with numerical results is shown in Fig. <ref> as the dot-dashed line. It appears that the D_1 and D_2 do not quite coincide with each other nor with Eq. (<ref>), which may be due to finite size effects.



It is instructive to draw an analogy between characteristics of the PBLMs and that of the Rosenzweig-Porter (RP) model from random matrix theory, see Ref. rosenzweig1960repulsion,Kravtzov2015RzPr and references therein, where the matrix elements are given by Gaussian random variable with zero mean and variance for diagonal and all off-diagonal matrix elements set ⟨ℋ_ii^2⟩=1 and ⟨ℋ_ij^2⟩∝ M^γ. Transition points between localized, delocalized and non-ergodic delocalized regimes as well as perturbative expressions for fractal dimensions Eq. (<ref>) are consistent in the two models. The dynamical correlator also shows similar behavior indicative of the formation of minibnads of non-ergodic eigenstates with the leading exponent 1-γ/2 in the scaling of the population transfer time with M coinciding in the two models. The prefactor (log M)^1/2 however is affected by the heavy tail of the distribution of the matrix elements and needs to be calculated analytically. It is  difficult to extract it accurately from the numerical simulations due to the finite size effects.





§  BORN APPROXIMATION FOR THE TRANSITION RATES



In this section we develop a simple picture relying on Fermi Golden Rule (FGR) to study the  rates of population transfer away from a given marked state to a set of other marked states inside the same miniband. 
Assume that the system is initially prepared at  a randomly chosen marked state |z_j⟩.  The probability amplitude  to remain in the initial state |z_j⟩ equals

    ψ(z_j,t)= ∑_βψ_β^2(z_j)  e^-i E_β t ,
 
where |ψ(t)⟩ evolves with the IB Hamiltonian  (<ref>) and |ψ_β⟩=E_β|ψ_β⟩. If  the eigenstates  dominantly coupled to the  marked state  |z_j⟩ are extended  then the amplitude ψ(z_j,t) will  undergo decay in time. 

Here we calculate ψ(z_j,t) using a simple effective Fano-Anderson model for the decay of a discrete state into a continuum <cit.>. This model captures the Born approximation for the ensemble of Hamiltonians introduced in Sec. <ref>. The model Hamiltonian   is obtained from the IB Hamiltonian  (<ref>) by  zeroing  out all off-diagonal matrix elements except those in the jth column and the jth row  connecting state |z_j⟩ to the rest of the marked states. The Hamiltonian  has the form

    = _j|z_j⟩⟨z_j|+ ∑_m≠ j (_m-iη)|z_m⟩⟨z_m|




    + ∑_m ≠ j_j m(|z_j⟩⟨z_m|+|z_m⟩⟨z_j|) ,

where the summation is over m ∈[1.. M],   m≠ j.  We consider the dynamics on a time scale when the population of the state |z_j⟩ decays  into the other states and introduce a  small imaginary part  -iη to  their energies. It is  assumed to be much bigger than the typical energy spacing,   η≫=W/M but smaller than the time scale on which the decay takes place.
 We introduce the parameterization similar to that in Sec. <ref> for the distribution of energies _j, 

    W=λ M^γ/2 ,

where λ is a (redundant) number of order of O(M^0).  
 

  








The amplitude ψ(t,z_j) has a well-known form <cit.>

    ψ(z_j,t)=∫_-∞^∞dz/π Σ_j^''(z)  exp(-i z t)/( z-Σ_j^'(z)-_j)^2+(Σ_j^''(z))^2 ,

where  we used a  short-hand notation 

    Σ_j(z)=_j(z)-i_j(z)

for real and imaginary parts  of self-energy of the marked state |z_j⟩

    Σ_j(z)=∑_m≠ j_jm^2/z-_m+iη ,

and we keep z real. Calculating the above integral to the leading order in   _j m (j≠ m) we get
 
    ψ(z_j, t)≃exp[-i(_j+Δ_j )t-1/2Γ_j t] ,

 where

    Δ_j≃_j(_j),  Γ_j ≃ 2_j(_j) .

 
The quantity Γ_j above is the total decay rate of the state |z_j⟩ which is twice the imaginary part of the self-energy  _j. The latter equals to the   "width" of the level _j due to the decay. Expressions (<ref>),(<ref>) correspond to a well-known Born approximation for the   self-energy _j. Using (<ref>) we get,

    _j=π∑_m ∈[1… M]/j_mj^2 δ(_j-_m,η) .

where we defined a function,

    δ(,η)≡1/πη/^2+η^2 .

The matrix elements _mj^2, see (<ref>),(<ref>), depend only on the Hamming distance d_mj.


 
 

The dominant contribution in to the sum (<ref>) comes from the transitions to the states  with |ϵ_j-ϵ_m| η. If the number of such states is large   the sum can be replaced by the integral corresponding to the approximation where the Lorentzian δ(ϵ_j-ϵ_m,η)≈⟨δ(ϵ_j-ϵ_m,η)⟩ is replaced with its average over realizations of ϵ_m. The average is independent of η≪ W which therefore drops out from the PDF of  the transition rate Γ and  the resulting  level width   =Γ/2.  
 This case corresponds to  the leading order Born approximation described  in the next Section <ref>.

  A more accurate treatment of δ(ϵ_j-ϵ_m,η)  as a random variable results  in the form of the PDF of Γ   (and ) being explicitly dependent on η. The physical meaning of η is the decay rate at the "children" sites ϵ_m, m≠ j which gives rise  to the  width  or the energy level _j at the parent site. In a large system the statistics of the decay rate for children and parents are expected to be the same. The crude approximation that captures this effect  is obtained by   substituting  η with typical value of . This corresponds to   self-consistent Born approximation described in  Sec. <ref>. It  gives rise to a more accurate expression for  PDF of  (and Γ) whose shape is rescaled compared to the leading order Born.
   Systematic analysis  is given by the cavity method described in Secs. <ref>, <ref>.
  
 



 §.§ Leading order Born approximation


We can break down the decay rate  Γ_j=2_j into a sum over different decay channels 

    _j=π∑_d=1^n V^2(d)(1-cos2ϕ(d)) ϱ_η^j(d) ,
 
where each term in the sum   corresponds to  the transition rate from  the initial state |z_j⟩  into  the subset of the marked states on a given Hamming distance d from  |z_j⟩ (see Fig. <ref>).  The factor ϱ_η^j(d)  in  (<ref>) is a spectral density of the marked states located at a distance d from the state 
|z_j⟩ within the window of energies η around _j 
 
    ϱ_η^j(d) = ∑_m≠ jδ(ϵ_j-ϵ_m,η)Δ(d-d_jm) ,

 where Δ(d) is a Kronecker delta and δ(,η) is defined in (<ref>). 



We denote  as M_j^(d) the number  of marked states that are separated by a Hamming distance d from the state |z_j⟩ (number of terms in the sum (<ref>))

    M_j^(d)=∑_m ≠ jΔ(d-d_jm) .

As discussed in Sec. <ref> the elements of the set   { M_j^(d)}_d=1^n  are   sampled from the multinomial distribution  with mean values 

    ⟨ M_j^(d)⟩= M p_d ,    p_d≃ 2^-nnd ,
where coefficient p_d  defined in (<ref>) is the probability  that a randomly chosen state is located on a Hamming distance d≠ 0 from |z_j⟩. The mean separation 
 between the adjacent energies _m in the sum (<ref>)  equals

    W/M_j^(d)∼  2^n/nd,   (M_j^(d)≥ 1)
    
     ,

where δ=W/M is mean spacing between the marked state energies. 
A substantial contribution to the sum in  (<ref>) comes from   the  terms 
corresponding to the marked states whose energy levels _j lie within the width η   from the energy _m, i.e., they  satisfy the resonant condition  |_j-_m|η as shown in Fig. <ref>.  

The contribution to a sum from each  resonance  is  ∼ 1/η and the
number of the resonances in a given decay channel  is ∼ M_j^(d)η/W (cf. Fig. <ref>). It is  shown in the Appendix <ref> that the dominant contribution  to the typical values of  _j (<ref>) comes from the values of d that correspond to ≫ 1. For them  the function δ(_j-_m,η) in Eq. (<ref>) changes  weakly between the adjacent values of _m,  and in the leading order Born approximation we estimate the  sum over m in (<ref>)  by  replacing it with an integral. 
Then  the spectral density can be estimated as 

    ϱ_η^j(d)≃ M_j^(d)  p(_j) ,

where we required  
 
    ≪η≪ W .

and p() is PDF of the marked state energies  with the width W (see  (<ref>)).



We plug (<ref>) into the expression (<ref>), obtaining  the following relation

    _j =π  p(_j)∑_d=1^n M_j^(d)V^2(d)(1-cos2ϕ(d)) ,

where the  sum  is dominated by values of d corresponding to  large values M_j^(d)≫ 1 (see Appendix <ref>). The steep exponential decrease  with d
of  the matrix element V^2(d) ∝ 1/nd (<ref>) is canceled by equally steep growth with d of the average number of states in the d channel ⟨ M_j^(d)⟩∝nd  (<ref>). As a result, the binomial factors cancels out and the average quantity ⟨ M_j^(d)⟩  V^2(d) changes only by  (n^-1) when d  changed by 1.


The term involving cos2ϕ(d) above  oscillates around 0 on the scale d∼ 1 (cf. Eq. (<ref>)). Therefore the contributions to the sum from the terms ∝⟨ M_j^(d)⟩cos2ϕ(d) average out.
 In what following we shall  neglect the cross-product of fluctuational and oscillatory parts  ( M_j^(d) -⟨ M_j^(d)⟩)cos2ϕ(d) and drop the second term  in the r.h.s of (<ref>) that contains cos2ϕ(d). 


Essentially, the above approximation corresponds to replacing the oscillatory part in the expression for the off-diagonal matrix elements  _ij≠ i=V(d_ij)√(2)sinϕ(d_ij) (<ref>) as follows:

    _i j→ V(d_ij) β_ij,  β_ij=±1,   i < j ,

where  β_ij are instances of a dichotomous random variable that takes values ±1 with probability 1/2. This approximate model of the ensemble of _i j will be also used in cavity method calculation in Sec. <ref>.

Using the expression (<ref>) for p()   and also Eqs. (<ref>), (<ref>),(<ref>),  we obtain the relation between the  PDFs  of the random variables  





    s_M,  =πV_ typ^2/W/M ,


    s_M=1/M∑_m=1^Mw_m .

Here  w_m are random variables independently sampled from the probability distribution g_∞(w) given in (<ref>). The level widths _j
  of  individual marked states for 1≤ j≤ M are  samples  of  the random variable  . 

In Eq. (<ref>)  we introduced  the characteristic value of the level width  . This equation relates 
  the PDF of   (or the decay rate Γ=2) to that of    and  M s_M.   
  We note that the  resulting  expression for the  level width  of a  marked state   formally corresponds to that given by FGR for the decay of the discrete level into the continuum <cit.>. The energies  of the marked states _m into which a given marked state |z_j⟩ decays form a miniband of the width _j. 
The  decay occurs simultaneously in many channels  corresponding to different Hamming distances between the initial marked state and the states of the  miniband.  


 The heavy-tailed PDF of the  random variable s_M is studied in details in Appendix <ref>. 
Using Generalized Central Limit Theorem (GCLT) for the  sums of a  large number of identical heavy-tailed random variables <cit.>
it can be represented in the form

    s_M  σ_M x+b_M

where   x obeys a  so-called Levy alpha-stable distribution  L^1,1_1(x)  <cit.> defined in the Appendix, Eq. (<ref>), and shown  in Fig. <ref>.
  Scaling factor and shift are

    σ_M=√(π/4log M) ,




    b_M≃σ_M^-1-2/πσ_Mlog(σ_M^-1) +2/π (1-γ_ Euler) σ_M ,

(γ_ Euler≃ 0.577 is the  Euler constant). They display very weak logarithmic dependence on M as compared with the main factor ∝ V_ typ^2/ in (<ref>). The width of the PDF  of s_M is shrunk   by a  factor (log M)^1/2≫1 and  the location of its maximum is increased
by  a factor (log M)^1/2≫1  compared to   L_1^1,1(x).






 The PDF of s_M has polynomial tail. Therefore  decay rates  of marked states Γ_j=2 _j can take a range values that are much bigger than their typical values 2 (<ref>), up to   M times bigger in the sample of the size M.  These atypically large decay rates   correspond to rare  clusters  of marked states that are located anomalously  close to each other. When clusters are formed by O(1) states the above  picture of the decay fails. 







§  NUMBER OF STATES IN A MINIBAND WITHIN BORN APPROXIMATION


 Using the  expression (<ref>) for the miniband  width  we can estimate the number of marked states Ω in a miniband corresponding to a given  state |z_j⟩. As before,  we divide the states  into the groups of the sizes , each corresponding  to the transitions away from |z_j⟩ with a fixed number of flipped bits d. The level width can  be written in the form  _j=∑_d=1^M_j,d where _j,d is  the  partial level-width  due to  the transitions  with flipping d bits.
 Then, using (<ref>)   and making use of the expression  (<ref>)  for the average values of M_j^(d) 
 we obtain 

    _j,d≃πV_ typ^2/δ 1/√(π n/2)M_j^(d)⟨ M_j^(d)⟩ .

The quantity  M_j^(d)/⟨ M_j^(d)⟩∼ 1 in (<ref>). This results in the interesting phenomenon   due to cancellation mentioned din the previous section: while   the typical number of marked states  in  a decay channel   
varies  very steeply with d,  typical values of  partial decay rates _j,d  in different channels do not.


The estimate for the typical number of states in the  miniband  at the distance d from |z_j⟩ is
 ∼_j,d/_j^d,

    ∼Ω  p_d,  Ω∼_j /∼(V_ typ/)^2 ,

where p_d=2^-nnd and Ω is the total number of states in the miniband.
 
One can also write the partial decay rate as Γ_j^(d)∼ V(d)Ω_d where the  product V(d)Ω_d does not depend on d (except from the prefactor). Of course, the analysis based on the decay rate does not apply for the  transition to the channels 
with  very few states. The condition
Ω_d≃ 1 leads to  Γ_j^(d)∼ V(d) for  d= corresponding to the typical Hamming  distance from |z_j⟩ to the nearest  marked state in a miniband 
where the condition V(d)≃_j^d  is satisfied (see Eq. (<ref>) in Appendix). 




The above  estimate gives the correct time scale ∼1/ V() over which the two states become hybridized. We note however that the total number of channels is  n-2=(n). As all Γ_j^(d) are nearly   the same, each channel contributes a small fraction (1/n) to the total rate. Therefore   V()∼Γ_j/n and  marked  state |z_j⟩  decays into the large number of marked states  within a miniband  before it has a chance to hybridize with the nearest one  at a distance .
 This property is markedly different from  the situation at finite dimension  <cit.>.
 
Using the scaling ansatz (<ref>) we estimate the mean separation between the energies   of marked states as

    =W/M=λ M^γ/2-1 .

Using the Eqs. (<ref>)  and (<ref>) we obtain the estimates for typical values of the decay rates and number of marked states in a miniband 

    Γ=2∼ M^1-γ/2,  Ω∼ M^2-γ .
 

We immediately observe that in the range of  γ>2 the number of marked states in a miniband  vanishes. It corresponds to a localized phase, consistent  with the fact that typical energy  spacing δ becomes  greater  than the typical tunneling matrix element  connecting the states. The number of states in a miniband   Ω cannot be greater than the total number of states M in the IB. The expression  above does not apply  for γ≤ 1.
This regime corresponds to ergodic phase.  

In the region 2>γ>1 the separation between adjacent eigenvalues of
 is of the same order as . The typical number of marked states in a miniband  Ω corresponds to   the typical number of non-ergodic delocalized eigenstates of  that form the miniband. 

    W≫Γ≫=W/M .

The number of states in a miniband scales as a fractional power of M less than one.  This is a hallmark of non-ergodic delocalized phase. 






§  CAVITY METHOD:  SUMMARY OF THE PREVIOUS RESULTS


 The cavity method has been actively used to study  Anderson Localization in  Levy matrices in the last several decades <cit.> starting from the seminal  work <cit.>.
  In the   present work we use cavity method to study the properties of minibands of delocalized non-ergodic states that were previously discovered in the studies of Rosenzweig-Porter  <cit.> and Regular Random Graph (RRG) <cit.> models.
  Initial studies suggested the existence of the mixed region with localized but non-ergodic states  <cit.>. However, recent numerical studies based on exact diagonalization using very large number of samples established that initially large crossover region between localized and extend states collapses in the limit of increasing matrix sizes <cit.>. Multifractal properties of eigenstates  in the localized phase and at criticality  were studied in  <cit.> using strong disorder perturbation theory.
 
Numerical solution of  cavity equations  to study localization transition in  Levi matrices with power-law distributions P(_ij^2)∝ 1/_ij^2(α+1) were obtained using population dynamics algorithm  <cit.>  utilizing the approach developed in  <cit.>. An alternative approach is based on the integral equation for the PDF of the diagonal elements of the resolvent <cit.>. It was obtained   in the limit where    imaginary part of the self-energy is vanishingly small <cit.>  (with the limit of infinite matrix size taken first). This allows one to derive  analytically the global density of states  <cit.> and the mobility edge E^*(α) which gives the 
 α-dependence of the energy   E^* separating extended and localized eigenvalues of  <cit.>. 

The cavity method proceeds as follows.  First,  we generate  a random M× M matrix  _ij (<ref>)  from the  ensemble  described in  Sec. <ref>. Then  we add  a new row (and a symmetric column) of independent numbers identically distributed as those in the old matrix  _ij. This is done by generating a random energy _0 from the distribution 1/Wp_A(/W); then  generating a random  bit-string z_0, computing  the array of Hamming distances d_j0 between z_0 and z_j and the corresponding matrix elements 
_j0=_0j for integer j∈ [1, M]. As a result we obtain a new  (M+1)× (M+1) matrix ^+1, where +1 emphasizes that it has one more row and one more column than . We will number elements of the new matrix by  indices running over the range [0,M] where the   index 0 corresponds to the added marked state |z_0⟩.  The cavity equations have the form    <cit.> 
    Σ_0^+1(z)=∑_m=1^M_0m^2G_m m(z) ,
   where 
    G_mm(z)=(z-_m- Σ_m)^-1 .
  It does not involve the non-diagonal matrix elements of the Green's  function G_mm^'(z) when statistical average ⟨_0m⟩=0. This is effectively our case as well (see Eq.(<ref>)).


The main assumption of cavity method is that in the limit  M→∞ the difference between the PDFs of Σ_0^+1(z) and Σ_0(z) disappears. This results in a self-consistent
 equations for the self-energy.  Following <cit.>  we add  small imaginary parts to the diagonal matrix elements _mm=_m-i η. It is a small “fictitious” quantity that is  still assumed to be much bigger than the marked state energy spacing 
  η≫ W/M. Results are  not expected to depend on the value of η provided its scaling with M is chosen appropriately, as will be discussed below.
  We separate the real and imaginary parts of the self-energy, Σ_m(z)=_m(z)-i _m(z) (cf. (<ref>)), obtaining
 


 
 
 
    _0     π∑_m=1^M_0m^2 δ(_m+η,z-_m-_m) ,
    _0     π∑_m=1^M_0m^2 δ(z-_m-_m,_m+η) .


where the function δ(x,y)≡1/πy/x^2+y^2 was already introduced in (<ref>).



The self-consistent  Eqs. (<ref>) were  derived  by Abou-Chacra, Anderson and Thouless <cit.> for matrices on Bethe lattices and by Bouchaud and Cizeau  for Levy matrices <cit.>.
The solution of these equation was only found in the case when they can be linearized in _m <cit.> giving the location of mobility edge E^*(α) as a function of the power α
in the tail of the PDF of the matrix elements P(_ij^2)∝ 1/_ij^2(α+1).
Here we will provide a full solution of the nonlinear equations.


We will solve the self-consistent equations (<ref>) under the assumption that pairs of variables (_m,_m) for each state m∈ [0,M] are taken from the same PDF 
(,;z) defined over the domain x∈(-∞,∞), y∈ [0,∞). In what following for brevity  we omit the explicit dependence on the parameter z. Following <cit.> we introduce the characteristic  function (k_1,k_2) of the PDF (,)

    (k_1,k_2)=∫_-∞^∞d∫_0^∞d (,)e^i k_1 +i k_2  ,

that satisfies the equation _η(k_1,k_2)=^M_η(k_1,k_2)



where








    (k_1,k_2) =
    
    
      e^i f k_1 δ(η+, z--)+ i f k_2  δ(z--,η+)
 
Here   f=_0m^2 and the average is performed with the joint PDF (,)1/W  p_A/WdfP(f).

The above relation  between (k_1,k_2) and (k_1,k_2) is actually an equation for  the PDF (,) because
 both  and  depend on  . 






§    SOLUTION OF CAVITY EQUATIONS  IN NON-ERGODIC DELOCALIZED PHASE




 §.§ Analysis of the imaginary part of self-energy


We note that the   exponent in the integrand of the above  expression for   depends on   and -z only via their combination +-z.
In the non-ergodic delocalized phase the typical width of the PDF of   is much more narrow than the width W of p() (<ref>).  We will also consider small values of |z|≪ W. Therefore in the first approximation we will neglect   and z compared to .  Then (k_1,k_2) depends only on the marginalized PDF

    ()=∫_-∞^∞d (,) .

Once this PDF is obtained, the PDF (,) can be analyzed from its characteristic
 function (0,k_2). Inverting it we obtain  the self-consistent equation for () in the limit M→∞

    ()=1/2π∫_-∞^∞dk e^M θ(k)-i k  .




    θ(k)=∫^∞_0df d dh   P(f) ()p_η+(h) (e^i k f h-1)

Here θ(k)=1-_η(0,k) and the  domain of integration for all variables is [0,∞). The function p_η+Y(h)  above is a conditional PDF of a random variable 
    h=δ(,η+Y)
 
with Y fixed and δ(x,y) given in (<ref>). The explicit form of  the PDF p_η+Y(h) is obtained  in  Sec. <ref> of the Appendix, Eqs. (<ref>),(<ref>).





To achieve further progress we use  the approximation (<ref>)  and drop oscillatory  factors in the off-diagonal matrix elements _0m. Then we have for the PDF
P(f)=g_∞(f/^2)/^2 (<ref>) and in what follows we will use the rescaled variable w=f/ ^2 for the squared matrix elements, in accordance with (<ref>). Instead of the variable h in (<ref>) we will use the re-scaled variable

    y=√(h(η+)) ,

 that  obeys the distribution 

    _η+(y)=2(η+)/W1/ y^2√(1-y^2)

(see details in Appendix <ref>, (<ref>)).
Then θ(k) takes the form

    θ(k)=∫_0^∞d() ϕ_+ηk ^2/+η .

Here ϕ_Y(u) is a characteristic function

    ϕ_Y(u)=∫_0^∞dx   g_Y(x)(e^i u x-1) .

of the PDF  g_Y(x) of the random variable x= w y^2 where  w obeys g_∞(w) and y obeys _Y(y) (<ref>).  Detailed study of g_Y(x) is given in  Appendix <ref>. The PDF  g_Y(x) depends on Y via the ratio Y/W and its plot is shown in  Fig. <ref>. It goes over into g_∞(y) for Y→∞.
 

We now make a key observation: in the limit of large x≫ 1 and for W≫ Y the following relations holds for the PDF  ϕ_Y(u)  and its characteristic function (see the corresponding Eqs. (<ref>) and (<ref>) in Appendix <ref>)

    g_Y(x)≃π Y/Wg_∞(x),  ϕ_Y(u)≃π Y/Wϕ_∞(u) .

The reason for this can be explained as follows.
For   large  deviations of   x=w Y^2/(^2+Y^2)   the conditional PDF p(|x) of the marked state energy   is  narrowly peaked in the range of values  ||∼ Y.
In contrast, typical  energy values are much bigger ∼ W.  This narrowing  of the conditional PDF p(|x)
gives rise to a small factor π Y/W in the  r.h.s. of (<ref>).


We observe that lim_k→∞θ(k)=0 and for M→∞ the integral  in (<ref>) is dominated by |k|≪ 1.
We make an  assumption (whose validity becomes obvious below) that for small enough k the integral in (<ref>) is dominated by values of  such that k^2/(+η)≪ 1. Therefore we will use in (<ref>)
the approximate expression for the characteristic  function ϕ_+η given by  Eq. (<ref>). We rescale  with the typical value of imaginary part of self-energy of marked states   (<ref>) obtained in FGR-based calculation in Sec. <ref>. Making a change of variables

    = s,  ()=1/ρ(s) ,

we rewrite the self-consistent equation (<ref>) in the limit x≫1, W≫ Y for the rescaled PDF ρ(s) in the following form:

    ρ(s)=1/2π∫_-∞^∞du e^-i u s+ Φ(u,Ω) ,




    Φ(u,Ω)=∫_-∞^∞dν ρν- β_η Ω ν ϕ_∞q/Ω ν ,


    β_η=η/ ,

  and

    Ω=π/δ=π/W/M^2 ,

 and  are  defined in (<ref>) and (<ref>), respectively. We observe that Ω corresponds to  the typical number of marked states in the mini-band  that we estimated in Sec. <ref> using  the Born approximation.

Assuming Ω≫1 (delocalized phase) we expand Ωνϕ_∞(q/(Ων)) in inverse powers of logΩ using asymptotic form of the characteristic function ϕ_∞(u) at small argument studied in Appendix <ref>,   Eqs. ((<ref>)),(<ref>). Truncating the expansion at terms ∼ (log M)^-1/2 we get 

    Ω  ϕ_∞(q/Ω)    ≃ -π |q|/2√(logΩ) ,
    Ω  ϕ_∞(q/Ω)   ≃  2q  (logΩ/π)^1/2+ q 1-𝒞-log |q|/(πlogΩ)^1/2 ,

where 𝒞≈ 0.577 is the Euler constant.
It is clear from comparing  individual terms in  Eq. (<ref>) with the exponential in  Eq.(<ref>) that q=(√(logΩ)). This justifies the order of truncation (see details in Appendix <ref>, Eq. <ref>). 

We make change of variables in the integral in (<ref>)

q=2 √(logΩ / π) t 
and obtain 

    ρ(s)=σ_Ω^-1 L_1^1,1((s-μ_Ω)/σ_Ω) ,

where quantity μ_Ω satisfies the   equation

    μ_Ω=b_Ω+2σ_Ω/π∫_-∞^∞ds ρ(s)log|s+β_η | .

Above   L^1,1_1(x) is Levy distribution    <cit.> defined in the Appendix, Eq. (<ref>), and shown  in Fig. <ref>.  Coefficients σ_Ω =√(π/(4logΩ))  and b_Ω≃ 1/σ_Ω are  given   in Eqs. (<ref>),(<ref>) where the   parameter M needs to be replaced by Ω. 

We plug the above expression for ρ(s) into (<ref>) and express μ_Ω in terms of a new variable x 

    μ_Ω≡ b_Ω-2σ_Ω/πlogσ_Ω^-1+σ_Ω x .

Then this variable   satisfies the following  equation

    x=2/π∫_-∞^∞ds  L_1^1,1(s)log |s+x+ζ_Ω| ,

that  involves a scale-free Levy distribution and  a single parameter ζ_Ω

    ζ_Ω=b_Ω/σ_Ω-2/πlog1/σ_Ω +1/σ_Ωη/

where  we used an explicit form of β_η (<ref>). We note that the  self-consistent equation for the   function ρ(s) is now reduced to the simple transcendental  equation (<ref>).

Using explicit form of   σ_Ω and b_Ω (<ref>),(<ref>) one can see that ζ_Ω is large compared to unity   in the delocalized phase,  ζ_Ω≃σ_Ω^-2 ∼ logΩ≫ 1. With  this property   the equation for x (<ref>) can be  solved by iteration using the asymptotic expansion of Levy distribution
at large arguments,  L_1^1,1(ν)≃ (2/π )ν^-2 (ν≫1). To the leading order

    x≃2/πlogζ_Ω+logζ_Ω/ζ_Ω .

Then using (<ref>) the expression for μ_Ω is

    μ_Ω≃1/σ_Ω+2σ_Ω/πlog1+η σ_Ω/+2σ_Ω(1-γ_ Euler)/π

where we neglected terms ∼σ^3_ΩlogΩ that are much smaller than the width σ_Ω of the distribution ρ(s)=σ_Ω^-1L_1^1,1((s-μ_Ω)/σ_Ω).

We note that the dependence of  μ_Ω (<ref>)  on the initial (fictitious) level broadening η disappears when the later is chosen to be much smaller  than the mini-band width <cit.>, W/M≪η≪σ_Ω.
 Using (<ref>),  (<ref>) the scaling behavior of η with M in the non-ergodic delocalized regime must satisfy the condition 

    η=M^κ,   |κ|< 1-γ/2,    γ∈(1,2) .



Finally  the expression for the distribution function of  the imaginary part of self-energy has the form

    ()=1/C L_1^1,1-/C ,


    =μ_Ω,   C=σ_Ω .

Here  is a   shift  of the distribution and C its scale parameter  (characteristic width). Also,

    μ_Ω   ≃1/σ_Ω+2σ_Ω(1-γ_ Euler)/π .
    σ_Ω   =√(π/4logΩ) .



Using the  scaling ansatz (<ref>) for the width W of the IB in terms of M ,  the typical number of states in a mini-band 
(number of resonances) equals, 

    Ω=π/λ^2 M^2-γ .

 
Using  the same scaling ansatz   (<ref>)  and  the expressions for  σ_Ω (<ref>) and μ_Ω (<ref>) we obtain,

    ≃2 π^1/2/λ  M^1-γ/2 (logΩ)^1/2 ,
    
    C   ≃π^3/2/2λ   M^1-γ/2   (logΩ)^-1/2 .








Here ∼n^1/22^-n/2 e^-n/(4 B_⊥^2)  is given in (<ref>).
The  shift  corresponds to  the typical value of .   One can see from the above that it is  logΩ∼log M≫ 1 times bigger than the distribution width. We note in passing that distribution of Σ^'' determines that of the miniband width Γ=2Σ^'' (<ref>).




  §.§.§ Comparison between the cavity method and leading-order Born approximation


It is instructive to compare the above distribution of  obtained using the cavity method with that obtained within the Born approximation (<ref>)-(<ref>).
In both cases the distribution of  is given by the appropriately rescaled and shifted Levy alpha-stable distribution L_1^1,1(x).
In both cases, the  scale parameter  C (characteristic width)  of the disribution  has the form C=σ_S with σ_S=√(π/(4log S)). In   the case of the Born approximation  S=M, corresponding to the total number of marked states, and in the case of cavity method S=Ω≪ M, corresponding to the (much smaller) number of states in the mini-band. Using (<ref>) we estimate

    σ_M/σ_Ω=√(2-γ) <1      (W=λ M^γ/2)  .

Therefore Born approximation underestimates the width  of the distribution of . The ratio (<ref>) is especially pronounced near the localization transition γ=2.  Value of σ_Ω^-1 shrinks to zero at the transition while that of 
σ_M  does not depend on the closeness to the transition point.

We note however that factors σ_Ω and σ_M depend on M only logarithmically.  At the same time, the leading order (power-law) dependence of the rescaling coefficient  on M is given by the factor ∝ M^1-γ/2,  and is identical in the cavity method and the Born approximation-based expressions. 

The situation is similar with the shift parameter   in the Levy distribution of  corresponding its  typical value,   ≃ /σ_S with S=M (Born approximation) and  S=Ω (cavity method). The  leading-order dependence of the shift on M is the same in both cases and is given by . In both cases the shift is greater than the rescaling coefficient by a factor ∼log M. However the Born approximation overestimates the shift by a factor  (2-γ)^-1/2 .




  §.§.§ Comparison between the cavity method and self-consistent  Born approximation




















The leading-order Born  approximation recovers    the typical shift  and the scale parameter  C  of the distribution of  with exponential accuracy in log M. However it gives an incorrect  dependence   of  the prefactor on log M in these coefficients.
The main approximation in Sec. <ref> was to assume that the sum in the expression for the spectral density ρ_η^j(d) 
(<ref>)  can be replaced by an integral. 
We revisit  the decay rate equation (<ref>) using the statistical ensemble (<ref>) 
 

 
 
 
    =^2/η∑_m=1^Mx_m,   x_m=w_m η^2/_m^2+η^2 .

Here in the l.h.s. we  omitted  the subscript in  _j and  made the rescaling V(d_jm)^2=^2 w_m. Random variables x_m  are sampled from the distribution g_η(x) given in (<ref>) and plotted in  Fig. <ref>). 
Using GCLT for the sum in (<ref>) one can  obtain the PDF of . The details are given in Appendix <ref> and here we provide the result,

    Σ̃_ typ^'' + x  C,

 
    Σ̃_ typ^''  = b_Ω_η  ,    C= σ_Ω_η  .

Here  x is a random variable that obeys Levy distribution L_1^1,1(x), coefficient 
σ_Ω  is given in (<ref>) and b_Ω is  given   in (<ref>) where one should replace M with the number of marked states in a mini-band of width η

    Ω_η= πη/δϵ.

Unlike the discussion  in the cavity method, the statistics of Σ^'' explicitly depends on η. 
We  make a self-consistent assumption and set η  equal to the characteristic  width of the miniband 

    η= ⟹ Ω_η=π/ .

 We conclude that the   typical number of states in a miniband Ω_η=Ω  given by the self-consistent Born approximation is the same as that  given by the cavity method,  Eqs. (<ref>). Therefore using (<ref>)  one can  see that the width   C  of the  distribution of   is also the same in both methods.
  The difference 
between the typical values of   in the two methods is  

    Σ_ typ^''-Σ̃_ typ^''=2/πC logσ_Ω^-1≪Σ_ typ^''

This error is much smaller than in the case  discussed in Sec. <ref> (cf. Eq. (<ref>)) where the self-consistent condition is not used.
However it  exceeds the distribution width C for sufficiently large   M≫1 because in the non-ergodic delocalized phase
 logσ^-1_Ω∼loglog M.
 



 §.§ Real part of self-energy


In this section  we will find the marginalized probability distribution of real parts  of self-energy 

    ()=∫_0^∞d (,) .

We  consider the first equation in (<ref>).  Following the arguments provided in Sec. <ref> we neglect the terms z-Σ^'_m in the r.h.s of the equation and   drop the oscillatory  factors in  _0m using the  probability distribution P(f)=g_∞(f/^2)/^2 (<ref>) instead. Then Eq. (<ref>) takes the form

    ∑_m=1^M r_m .
 
Here r_m are instances of a random variable R such that

    r= f /^2+ ()^2 ,

where  , f, Σ^'' are random variables independently sampled from the distributions p(), P(f) and (), respectively. Using GCLT, in the asymptotic limit of M→∞ the sum in (<ref>) is determined by the tail of the probability distribution  of r at |r| →∞. This analysis is very similar to the one already discussed in Sec. <ref>,<ref> and in Appendix <ref>. Here we omit details of the calculations and simply provide the result.
The tail of the PDF of r in the limit |r|→∞ has the form

    ρ=r/2/(π M),   PDF(ρ)≃1/ρ^2√(log(ρ)/π) .

(ρ≫ 1). The distribution function () of the sum in (<ref>) is the  Cauchy distribution

    ()=1/π/()^2+()^2,  =/σ_M .

Here the expression for σ_M∼ 1/√(log M) is given in (<ref>). Cauchy distribution has the form very similar to the stable distribution L_1^1,1(x) that describes the fluctuations of the  (<ref>) up to the shift and rescaling coefficients. Both distributions are displayed in Fig. <ref>. The tail of the Cauchy distribution differs  from that of L_1^1,1(x) by a factor of 2. Unlike that of  the  distribution of 
  is symmetric for impurity states with energies near the center of the band. 
The typical value of   is greater than that of  by a constant factor

    / =1/√(2-γ)     (W=λ M^γ/2) .

The width  of the distribution of  is the same as its typical value while the width  C of the distribution of  is   smaller by a factor ∼  1/log M (cf. Eqs. (<ref>),(<ref>)).  These relations between the distributions of  and  have implications for the complexity of the population transfer as discussed below. We also note that the real and imaginary parts of self-energy of  a given marked state are correlated with each other because according to Eqs. (<ref>),(<ref>) the values of  _j and _j depend on the same set of parameters (_jm, _m, etc). In this work we will not study their correlations.




 §.§  Dynamic correlations


For states close to the center of the band of marked states the typical value of the mini-band width can be connected to the average of the dynamical correlator, with the delta function regularized by a finite scale η, Σ”_ typ≫η≫δϵ, δ(x) →δ_η (x)≡1πηx^2+η^2, 

    1/Σ”_ typ+η=1/π∫ _-∞^∞dωη/η ^2+ω ^2   p (ω),

which can be inverted to obtain,

    p ( ω) ≈[
     [                1/Σ”_ typ,                    ω≤ω_Th; 1/Σ”_ typ(ω _Th/| ω| )^2,                    ω>ω_Th ].

where we introduced the Thouless energy,

    ω _Th=1/2πΣ^''_ typ.

The typical value of the mini-band width was obtained in Eq. (<ref>). From the comparison of the respective Fig. <ref> we conclude that the scaling of the typical population transfer time 1/ω_ Th and the scaling of the value of the dynamical correlator K(ω) are consistent in numerical and analytical calculations, subject only to a small correction in the scaling exponent ε=0.05.





§  COMPLEXITY OF THE POPULATION TRANSFER PROTOCOL


After the system is prepared at a given marked state |z_j⟩ at t=0 the probability for the population to be transferred to other marked states is 1-ψ^2(z_j, t). At the initial stage 
the survival probability decays exponentially (<ref>) with the mean decay time 1/Γ_j=1/(2_j).


The  initial marked state decays into the eigenstates |ψ_β⟩ of the IB Hamiltonian  with typical energies E_β inside the narrow interval  corresponding to the miniband associated with |z_j⟩. It has a width  _j and is centered around _jj=_j.  Typical classical energies  of the bit-strings measured at the end of PT protocol will obey the probability distribution (-_j-_j) with  given in (<ref>). The  success of PT protocol  is  to find a bit-string distinct from z_j at a time t with energy inside that window  Δ_ cl around _j.  The expected time to succeed in PT equals

    t_ PT^j= 1/2_j p_Δ,   p_Δ=∫_0^Δ_ cl-_j-Δ_ cl/2d .

Here p_Δ is the probability of detecting a bit-string inside the target window  Δ_ cl under the condition that initial state has decayed.
Let us assume that the PT window  is as wide as the typical miniband  width, Δ_ cl=.  In this case p_ mb differs from 1 only by a constant factor that does not depend on M (cf. (<ref>)). Therefore we will detect the bit-string inside the PT window with finite  probability as long as we waited long  enough for the  transition away from the initial marked state to occur. Because the initial state  |z_j⟩ is picked at random we can estimate typical  time  to success  of PT   ∼ 1/  corresponding to the inverse typical  width of the miniband.   All of the states in a miniband are populated at (roughly) the same time  because    transition rate to a subset of states on a distance d away from |z_j⟩ depends on d very weakly (see Eq. (<ref>) and related discussion in Sec. <ref>)).
 
 From a computational perspective it is of interest to characterize the PT  by the relation between the  the typical success time   of PT   and the number of states Ω over which the population is spread  during PT

    ∼1/ √(ΩlogΩ)∼(2^n/n ΩlogΩ)^1/2 e^2θ n .

where we set Δ_ cl∼ (see discussion above).  We note that the time t_ G for the Grover algorithm for unstructured quantum search to find Ω items in a database of the size 2^n is t_ G∼ (2^n/Ω)^1/2. 
PT time  scales worse than Grover time t_ G by an additional exponential factor e^2θ n≃ e^n/2B_⊥^2   (<ref>). The  scaling exponent 2 θ  can be made arbitrarily small at large transverse  fields 1≪=(n^0).
 
One can expect that the distributions of eigenvalues and eigenvectors inside the mini-band  are very similar to those    in the ergodic case, albeit with the appropriately rescaled effective dimension Ω of the Hilbert  space <cit.>. For example,  the energy spectrum of the  mini-bands in the non-ergodic delocalized phase of Rosenzweig-Porter (RP) model corresponds to  the Gaussian Orthogonal ensemble.
There, according to  the semicircle law <cit.>,  the typical spectral width of the mini-band  (∼ 1/)  is proportional to the square root of the number of states Ω in it. 
Therefore the Grover scaling  (<ref>) for PT is consistent with semicircle law in the Gaussian random matrix models  that allow for non-ergodic delocalized phase such as RP model.


 However in the case of Levy matrices the distribution of eigenvalues has polynomial tails <cit.>,  their spectrum is not bounded and semi-circle law does not apply. As mentioned above, this leads to a broad distribution of PT rates.  There exist statistically significant clusters of states of a relatively small size that will be populated faster than typical case because the corresponding classical bit-strings  are located closer  to each  in Hamming distance than the  typical inter-state separation. At first glance,  this tendency is counter to the Grover scaling (<ref>). We note however that fluctuations of  and  are correlated with each other. Faster  decay of a marked state will also correspond to bigger self-energy shift which will reduce the likelihood of finding a marked state with its energy inside the target window Δ_ cl∼.
 
However the Grover  scaling still survives in a typical case corresponding to PT away from a randomly selected bit-string. For Levy matrices  <cit.>  it reflects  the fact that  the typical width  of the curve of the global density of states along the energy axis must scale as a square root of the corresponding typical number of states (area under the curve).





§  COMPARISON WITH THE ANALOGUE  GROVER SEARCH




 §.§  Grover search starting from a fully symmetrized  state


So far we have studied the PT protocol with the Hamiltonian (<ref>) H=H_D+H_ cl that  starts from a  given marked state of an IB model  H_ cl (<ref>)  and aims at finding a different  marked state inside a given window of energies using a transverse field Hamiltonian H_D=-∑_m=1^nσ_m^x (<ref>) as a driver. 

We  consider here a different protocol inspired by the Hamiltonian version of Grover algorithm proposed in <cit.>. The new protocol 
finds  marked states in the IB model H_ cl starting from the  ground state of H_D   which is a
fully symmetric state |S⟩=2^-n/2∑_j=1^n|z⟩  in a computational basis. This protocol can be implemented   by adjusting  the value of transverse field  ≈ 1 so that  the ground state energy of the driver is set near the center of the IB. Then we can replace the full driver with the projector on its ground state, H_D→ -n|S⟩⟨S|. 
The quantum evolution is guided by the  Hamiltonian: 

    H_ G=-n|S⟩⟨S|+∑_j=1^M(z_j)|z_j⟩⟨z_j| .

With the   initial condition  |ψ(0)⟩=|S⟩. In the case where all impurity energies are equal to each other, {(z_j)=-n}_j=1^M, and =1 the Hamiltonian H_ G  is a generalization of the analog version of Grover search <cit.> for the case of  M target states. The system performs Rabi oscillations between the initial state |S⟩ and  the state which is an equal superposition of all marked (solution)  states. Time to solution is the half-period of the oscillations, the "Grover time" t_ G

    t_ G=π/2n√(2^n/M) .

Hamiltonian versions of Grover search with transverse field driver whose ground state were tuned at resonance with that of the solution state were considered in <cit.>.

Robustness of the Grover algorithms to phase noise was considered previously in the case of a single marked  state <cit.>. Here we investigate the role of 
systematic phase errors in quantum oracle for the case of multiple  solutions by assuming that  marked state energies   take distinct values (z_j)=-n+_j
randomly   distributed over some  narrow  range  W. We will also investigate the systematic error in the Grover diffusion operator<cit.>. In the Hamiltonian formulation <cit.>
this corresponds to the deviation from unity of the parameter   that controls the weight of the driver in (<ref>). We will define

    =1-_0/n ,

where  _0 is the driver error. 

We denote the computational basis states as |j⟩≡|z_j⟩ with j∈[1,N],  N=2^n and assume that marked states correspond to the range  j∈[1,M]. We also introduce the state 
 
|0⟩=1/√(N-M)∑_j=M+1^N|j⟩

that  is  orthogonal  to all  the marked states. The subset of  basis vectors ={|j⟩}_j=0^M spans the M+1 dimensional subspace  with the remaining set _⊥ of  basis vectors spanning the  orthogonal N-M-1 dimensional subspace. One can show that  H_G does not have matrix elements that couple   with _⊥. 




Assuming that N≫ M one can  consider the  decay of the  state |0⟩ instead of the state |S⟩. 
We use (<ref>) and omit constant terms and small corrections (M/N) in H_G.  The  non-zero matrix elements  H_G^ij=⟨i|H_G|j⟩ 
in this subspace  have the form 




    H_G^jj=_j,    H_G^j0= -(1-δ_j0)V,   V=n 2^-n/2 ,

where  j∈[0,M] and H_G^j0=H_G^0j.
On a time scale t≪1/=M/W much smaller than the inverse spacing of the  energies _j
the quantum evolution with  initial condition |ψ(0)⟩=|0⟩ corresponds to the decay of the discrete state with energy _0 into the  continuum <cit.> with the finite spectral width W <cit.>. It is a similar problem to that  discussed in the Sec. <ref>.







  



  §.§.§ Sensitivity to  systematic oracle phase error


We first consider the case of  relatively large oracle errors (wide energy band W) 

    V√(M)≪ W ≪ V M ,

and modest driver errors 

    _0=n(1-)  W .

In this case, following the results of the Sec. <ref> on the  solution of the Fano-Andreson model  <cit.> we obtain an exponential decay of the initial amplitude
(cf. (<ref>))

    ψ_0(t) ≃exp[ -Σ_0^'' t-i_0 t-i Σ^'_0(_0+i 0^+)t ] .

where Σ_0(z)=Σ_0^'(z)+i Σ_0^''(z)  is a self-energy
and 

    Σ_0(z)=V^2∑_m=1^M1/z-_m,  Σ_0^''≡1/2Γ_0=π V^2/W/M .

The state |0⟩  undergoes an exponential decay  with the  rate Γ_0=2Σ_0^''. After the characteristic  time ∼ 1/Γ_0  the population is transferred into a  subset of the marked states with energies inside the window |_j-_0|≃Σ_0^''≪ W.

 The number of marked states (solutions)  to which the population is transferred is Ω∼Σ^''_0/.
 The relation between   and  Ω is
 
    ∼1/V√(Ω),  Ω∼(V/W/M)^2 ,

 the same  as in the Grover algorithm (<ref>). It also recovers the scaling with Ω and n, up to a factor exp(-n/(2^2)), for the  time of PT  considered in the rest of the paper that uses  transverse field as a driver  and starts from any  marked state instead of a fully-symmetric state. 
 
 
 To characterize the effect of oracle errors we introduce the scaling  ansatz  for the marked states bandwidth W∼ 2^-n/2 M^γ/2 similar to that in   (<ref>).  We observe that the number Ω of solution states populated over the time   cannot be greater than M by construction.  For  W V√(M) (or γ<1) the value of Ω≃  M and the scaling of the transfer time  with M is the same as t_ G in the ideal  Grover  algorithm (<ref>). In the region given by  (<ref>) (or 2>γ>1)  the algorithm performance is degraded because Ω≪ M. For    W≫ V M (or γ>2) 
 the algorithm fails to find even one solution. 
 



  §.§.§ Sensitivity to  the systematic driver error



We now consider the sensitivity of the algorithm to an error in  the  weight of the driver Hamiltonian, i.e., to the nonzero value of the parameter e_0=n(1-) (<ref>).  We assume that 
_0≫ W while the spread of the marked state energies

  the condition (<ref>), so that absent driver errors, PT time  would follow a Grover-like scaling law (<ref>).


In this case the state |0⟩ is coupled non-resonantly to a continuum with narrow bandwidth.  The   expression for the population transfer  to the marked states can be obtained  from the time-dependent perturbation theory in the parameter _0/W

    ∑_m=1^M|ψ_m(t)|^2=2 M V^2/_0^21-cos(_0 t)sin( Wt/2)/W t/2) .
    
    Maximum transfer occurs at the time t_0=π/_0 with the total transferred probability p_0=4M V^2/_0^2. Typical time ≃ t_0/p_0 to achieve the successful population transfer to marked states involves repeating the experiment 1/p_0 times 
    
        =1/Γ_0 π^2_0/W ,
    
    where Γ_0 is given in (<ref>) and the first multiple in r.h.s gives the typical transfer  time in the absence of driver errors. The later leads to an increase of the transfer time by a large factor _0/W.
    
    
    For  the maximum possible bandwidth  W when nearly all states are populated, W∼Γ_0∼V√(M),   the time of population transfer (<ref>) is 
    
    
    
    
        ∼ t_ G(t_ G_0 )     (_0 ≫ t_ G^-1∼ V√(M)) .
    
    As expected, when  the driver error exceeds inverse  Grover time 1/t_ G the performance of analogue Grover  algorithms  (<ref>)  degrades  relative  to t_ G. 
     This is a direct consequence of the fact that the quantum evolution  begins from fully symmetric state which is a ground state  of the driver Hamiltonian whose energy is tuned at resonance with the marked states. In this case  the transverse field Hamiltonian driver effectively corresponds to  the projector (<ref>). Because the ground state is not degenerate, the resonance region is exponentially narrow (∼ 2^-n/2√(M)). This  results in the   exponential  sensitivity of the Grover algorithm performance  to the value  of driver weight. This critical behavior was studied in the work on quantum spatial search <cit.> for the case of one marked state.
    
    In contrast, in the PT protocol considered earlier in the paper there was no need to fine-tune  the value of  other than making  it  large, ≫1. 
    This happened  because   the effective coupling between the marked states described by the down-folded Hamiltonian  (<ref>) was not due to any one particular eigenstate of the driver (such as the state |S⟩ for the Grover case). Instead this coupling was 
    formed due to an exponentially large (in n)  number of  non-resonant,  virtual  transitions between the marked states and highly  exited states  of the transverse field Hamiltonian H_D.  This resulted in a significant  improvement in robustness for the proposed PT relative to the analogue Grover  algorithm.  
    
    
    
    
     §.§  Grover search starting from a marked state
    
    
    We now consider an implementation of the analogue Grover search that starts from the marked state similar to the PT protocol considered in previous Sections. The transition amplitude U_ij(t)=⟨i|exp(-i H_ G t)|j⟩ between the two marked states  can be written in the form
    
    
    
    
        U_ji(t)=∑_λ e^-iλ tψ_λ(i)ψ_λ(j) .
    
    Here  ψ_λ(j)=⟨j|ψ_λ⟩ are amplitudes of the eigenstates of H_ G in  the M+1 dimensional subspace and λ are
    the corresponding eigenvalues  that obey the equation 
    
        λ=_0+∑_j=1^MV^2/λ-e_j,  ψ_λ(j)=V/λ-_j1/√(Z_λ) .
    
    Here
    
        Z(λ)=1+∑_m=0^MV^2/(λ-_m)^2 .
    
    Instead of providing a detailed analysis of the above solution we  provide  an order of magnitude estimate to extract the relevant scaling behavior. 
    We again assume that   the spread of the marked state energies, W=t_ G^-1=(V √(M)) corresponds to the inverse of the Grover time  t_ G needed  to find any one of the solutions with equal probability.  The typical separation between the adjacent vales of _j is  =W/M∼ V/√(M). 
    
    It follows from (<ref>) that in the  ordered array obtained by  combining  together the sets of  energies {e_j}_j=0^M  and eigenvalues {λ_m}_m=0^M their values  appear alternatively and sequentially, e.g., _j-1<λ_j<e_j<λ_j+1. The typical separation between the adjacent elements   in the array is |λ_j-_j|∼. We observe that for a given value of λ the  sum in the expression for Z(λ) (<ref>) is dominated by  the small, (1),  number of  terms  with  |_m-λ |∼, each term of the order of M.  Indeed, there are  (M) remaining terms   corresponding to   |_m-λ |∼ W. The  magnitude of those terms is V^2/W^2∼ 1/M and their aggregated contribution to the sum is (1). Therefore we can estimate  Z(λ)=(M) and for the amplitudes we have
    
        ψ_λ(m) ∼V/λ-_m1/√(M),    m=i,j .
    
    For a given initial state |i⟩ at time t we pick the final state |j⟩ within the energy window _j-_j∼Δ = 1/t around _i. The sum in the expression (<ref>) for the transition amplitude U_ji(t) is dominated by the number of terms Ω= Δ/∼Δ√(M)/V corresponding to the  eigenvalues λ inside  the same  window of energies. For those terms  λ-_i, _j-λ∼Δ giving the estimate for the amplitudes ψ_λ(i), ψ_λ(j)∼ 1/Ω (cf. (<ref>)). The magnitude of the sum in (<ref>) can be estimated as |U_ij(t)|∼Ω|ψ_λ(i)ψ_λ(j)|∼1/Ω. On the other hand, because ordered values of λ and _m alternate in sequence the probability |U_ij(t)|^2 is distributed over Ω marked states and |U_ij(t)|∼Ω^-1/2. By equating the  above two estimates for |U_ij(t)| we immediately obtain Ω∼ 1 and therefore
    
        Δ =1/t∼∼V/√(M),   (Ω∼ 1).
    
    
    In the case when there are only a few marked states (M∼ 1 and W∼ V)  the  
    probability is initially localized on a given marked state |i⟩ and then it spreads  over to others states  separated in energy by V during the time  t_ G∼ 1/V∼ 2^n/2. In this case the algorithm 
     time scales with n identically to that of the analogue Grover search that starts at the fully symmetric state |S⟩. Similar performance is achieved by the PT protocol using transverse field ≫1 and discussed
     in previous sections. 
    
    The difference from analogue Grover search starting at |S⟩ from the above PT protocol using a transverse field becomes  dramatic for large number of marked states M≫1. Both analogue Grover search and the PT protocol  benefit from the increase in M: the algorithmic time shrinks ∝ 1/√(M) and the number of marked (solution) states Ω  in the number of states in the final superposition increases  with M.
    
    In contrast, the quantum search with H_ G starting form the marked state |i⟩ does not create massive superpositions of solution states when M increases. Instead it involves  a very few others states that are adjacent in energy, |_j-_i|∼ V/√(M).  The time of the algorithm increases with M (<ref>). This happens because unlike the Hamiltonian H with a transverse field (<ref>), the Hamiltonian H_ G is integrable. The wave-function remains localized near the initial marked state.
    
    
    
    
    § OUTLOOK
    
    We believe that the  approach developed in this paper can be generalized beyond the IB model to  study PT  in other types of transverse field quantum spin-glass problems (<ref>).
     In the typical problem of interest, the classical energy landscape  contains a large number of deep local minima separated  by  Hamming distances d=(n) and an exponential number of them, M=2^μ n, μ=(1), are located  within  some narrow strip  Δ_ cl  of  classical energies   chosen at the tail of the distribution P_ cl() that decays exponentially with the distance from its maximum.   Such situation occurs, e.g.,  in the Random Energy Model <cit.>,  p-spin model with p>2<cit.>,   and other problems under certain choices of the ensemble parameters. 
    
    The objective as before is to estimate the statistics and the  scaling of the miniband sizes with n in the  energy strip  Δ_ cl  for large enough values of n, where finite size effects are small, typically n∼ 100 or more. Direct numerical solution of the Schrödinger  equation is no longer computationally tractable at these sizes. 
    
    
     
    Essentially,  our  approach is based on the existence of  two distinct energy scales. The first scale is the  typical change in classical energy  corresponding to one bit flip: _ flip. The second scale is the typical width of  non-ergodic minibands Γ < Δ_ cl, which decreases   exponentially with n.
    The tunneling transitions between the states inside the  miniband  require  a large number of spin flips, and therefore 
    _ flip≫Γ. Starting from the initial state |z_i⟩ inside the strip of energies Δ_ cl,  the quantum evolution   is confined within the corresponding miniband. This  evolution can be described by an effective down-folded Hamiltonian _ij
    defined over a subset of computational basis states whose classical energies lie within the energy strip Δ_ cl.
    
    
     The effective Hamiltonian  arises in the  the self-energy formalism <cit.>. The   off-diagonal matrix elements _ij can be understood in terms of the locator expansion  <cit.>, i.e, by expanding the   matrix elements of the resolvent  of H (<ref>)   in powers of the transverse field . Each matrix element 
     _ij  corresponds to the sum over all consecutive  spin-flip processes  that  begin in the state |z_i⟩, proceed through virtual states separated by energies at least _ flip from the corresponding miniband, and return  to the miniband only at the last step, at the state 
    |z_j⟩.    
      We emphasize that _ij takes into account the trajectories containing loops, i.e., multiple visits of the high-energy virtual states.
    
    
     
      
    
    
    The off-diagonal  elements  _ij are exponentially small in n  since the  Hamming distances between the states are extensive, d_ij= O(n). They are determined by the elementary spin-flip process  of a  very high order  in . The distribution of  _ij can be estimated either  analytically via the WKB approach or  numerically via the forward scattering approximation <cit.>. 
    
    
     At the same time the renormalization of the diagonal matrix elements _ii-(z_i)  appears already in the second order in . This fact is not important for the IB model where the renormalization is exactly the same for all marked states z_i and thus leads to a uniform shift of their energies. Therefore if _ii and _jj are separated by a very small energy then  so are   (z_i) and (z_j).
     
    Unfortunately generally  this is not the case.  Repeating PT experiments  starting from the same state  |z_i⟩ will produce  a set of  states  |z_j⟩ at the output.
    The dispersion of the corresponding values of _jj-_ii  is very small – of the order of the miniband width Γ_i.
    However, the dispersion of the classical energies _j is much  bigger due to renormalization of the diagonal matrix elements _ii-_i by transverse field. The   values of _j observed in PT will be distributed with the  width (n^1/2) around some mean value shifted by (n)  from the center of the miniband (see Fig. <ref>). Therefore  one must consider the PT window that is much wider than the miniband width, Δ_ cl=(n^1/2). 
    
    
    Suppose that  a classical algorithm proceeds via a continuous-time random walk between the states inside the window Δ_ cl with the transition rate matrix Q_ij and a typical matrix element Q_ typ.  Given the  PDF of classical energies  P_ cl() the separation between the neighboring energy values is =1/(2^n P_ cl(_j)) and the number of states with the classical  energies inside the strip Δ_ cl is 2^n P_ cl(_j) Δ_ cl.  Then the quantum () and classical (t_ cl) search times starting from |z_i⟩ are  
    
        t_ PT∼V_ typ^-2/2^n P_ cl(_i),   t_ cl∼Q_ typ^-1/2^n P_ cl(_i) Δ_ cl ,
    
    where V_ typ is the typical value of the off-diagonal matrix elements. Here we used the estimate (<ref>) for the quantum transition rate that also determines the width of the mini-band Γ_i∼1/.  Assuming  that both the quantum and classical searches are   unstructured we get  V_ typ^2 Q_ typ∼ 2^-n. 
    It follows from Eq. (<ref>) that 
    the quantum speed-up observed in IB model
    is lost because Δ_ cl scales down at most polynomially with n. 
    
    
     However the  PT computational primitive  can be modified in a way  that the quantum speed up is restored. Given an initial bit-string z_i the task is to find a bit-string z_j≠ z_i such that |_ii-_jj| ≤Δ_ cl.
      Quantum PT protocol to solve this problem is to start from the state |z_i⟩,  evolve the system under the Hamiltonian H(<ref>), measure in computational basis,  and  check if the outcome z_j is the solution.
    
    Because the dispersion of the renormalized energy values  _jj is small for the states within the   same miniband,  the width of the PT window can be chosen very narrow  Δ_ clΓ_i∼ 1/. Consequently,   the quantum speed up in Eq. (<ref>) is  recovered. 
     
     It is crucial that we can  check  the  correctness of the  solution with a classical computer: given the bit-string z_j  a classical algorithm should  determine the value of _jj  with the accuracy  (Γ). This could be done by  computing numerically the locator expansion of the resolvent of H (<ref>)  in powers of . An alternative approach is to use  Stochastic Series Expansion (SSE) algorithms <cit.> to compute the  contribution to the  partition function Tr exp(-β H)  from the  the valley of (z) containing the bit-string z_j by expanding it in powers of inverse temperature β and transverse field <cit.>.[This can be done 
     if the   state |z_j⟩  is sufficiently close in energy to the corresponding deep local minimum of (z).
     
     ]  We note that the level of precision required  does not involve the  hybridization of the state |z_j⟩ with the states  from other deep minima with comparable energies (cf. Fig. <ref>). This    limits the Hamming distance between  the  state |z_j⟩ and the computational states involved in the SSE computation.  The relevant number of terms of  the aforementioned   expansions as well as the computational time scale with n  as
      1/Γ∼, i.e., inversely proportional to the required precision.
      
     
      
    
    
    We note that for the bit-string z_i  close to  a local minimum of (z),   Quantum Monte Carlo (QMC) could simulate a tunneling transition away from z_i with the  rate Q_ij that scales   as  Q_ij∼_ij^2<cit.>. Therefore, from Eq.(<ref>) we have
    
    t_ cl^ QMC∝ t_ PT^2,
    
    which corresponds to a quadratic speedup of  the  quantum PT algorithm over QMC for this particular sampling task.  In general,  the states within   the narrow energy belt Δ_ cl are far away from the local minima  and the QMC methods cannot be directly applied.  In summary,  we argue that coherent multi-qubit tunneling can be a powerful computational resource for solving PT problems with complicated  energy landscapes.
    
    
    
    § CONCLUSION
    
    
    
    We analyze
    the computational role of coherent multiqubit tunneling that gives rise to  
    
    bands of nonergodic delocalized quantum  states as a coherent pathway
    for population transfer (PT) between computational states with close 
    energies.  In this regime PT cannot be efficiently simulated by QMC.  
    
    We consider optimization problems with an energy function (z) defined over the set of 2^nn-bit-strings z. 
    We define a computational primitive  with the objective to find   bit-strings  z_j≠ z_i  inside some narrow energy window Δ_ around the energy of the initial bit-string z_i.
    The problem is hard for sufficiently low   starting   energy (z_i) in the region  proliferated by  deep local minima   that are separated by  large Hamming distances.   
    
    
    
    
    We propose to solve this problem  using the  following quantum population transfer (PT) protocol: prepare the system in a computational state |z_j⟩ with classical energy (z_j), then evolve it with the transverse-field quantum spin Hamiltonian. Classical energies (z) are encoded in the problem Hamiltonian diagonal in the basis of states |z⟩ similar to quantum annealing (QA) approaches <cit.>. A key difference  from  QA or analogue quantum search Hamiltonians   <cit.> is that the transverse field is kept constant throughout  the algorithm and  is not fine-tuned to any particular value.   At the final moment of PT  we
    projectively measure in the computational basis and check if the
    outcome z is a “solution”, i.e.,  z  z_j, and the energy
    (z) is inside the window Δ_.
    
    
      In this paper we  analyzed PT dynamics in  Impurity Band (IB) model with a “bimodal” energy function: (z)=0  for all states  except for  M“marked”  states  |z_j⟩ picked at random with energies forming a narrow band of the width W separated by a large gap (n) from the rest of the states.  This landscape is  similar to that in analogue Grover   search   <cit.> with multiple target states and a distribution of oracle values for the targets. The best known classical algorithm for finding another marked state has cost O(2^n/M). 
    
    The transverse field gives rise to tunneling between a pair of marked states  corresponding to a  sum over a large number of virtual transitions connecting the two marked states via the states with (z)=0.  As a result the PT dynamics  is described by the down-folded
    M× M Hamiltonian  that is dense in the space of the marked states |z_j⟩. Its  off-diagonal matrix elements  _ij=V(d_ij)cosϕ(d_ij) depend only on the Hamming distance d and are obtained using  WKB method. The distribution of matrix elements _ij  has a heavy tail decaying as a cubic power for V(d)≫ V_ typ. This is a remarkable result of the competition between the very steep decay of the off-diagonal tunneling matrix element with the Hamming distance d, and the steep increase in the number of marked states M_d∝nd at distance d. We emphasize that such polynomial tail in the distribution of matrix elements is only possible either in infinite dimension or in presence of long-range interactions (e.g, dipolar glass). 
    
    The dispersion of the diagonal elements  _jj=(z_j) is expected to be large, W∼ V_ typ M^γ/2≫ V_ typ with γ∈[1, 2]. Therefore we call  _ij a Preferred Basis Levi matrix (PBLM), a generalization of the Levi matrix from the random matrix theory. We demonstrate two localization transitions in the PBLM ensemble whose locations are determined by the strong hierarchy of elements of the PBLM  _ij. In the range 1<γ<2 there  exist  minibands of non-ergodic delocalized eigenstates of . Their  width is proportional to 1/≪ W.   Each miniband associated with a support set  over the marked states.   If  γ>2 then  W exceeds the largest matrix element of  _ij and the support set is empty, all  eigenstates  are localized. If γ<1  then  W is smaller than the typical largest element in a  row of _ij and the support set extends to all marked states – all eigenstates   are “ergodic”.
    
    
    We find the distribution of the miniband width Γ=1/ analytically by solving the non-linear cavity equations for an ensemble of PBLMs. Unlike previous analyses focused on linearized cavity equations near the Anderson transition, we find the solution of the fully non-linear cavity equations in the non-ergodic delocalized phase. 
    
    The distribution of miniband widths Γ  obeys alpha-stable Levi law with tail index 1. The typical value of Γ and its  characteristic variance  exceeds the typical matrix element of  by  a factor Ω^1/2 where Ω=(M V_ typ/W)^2 is a size of the support set in a typical miniband.
    
    We demonstrate that quantum PT finds another state within a target window of energies Ω in time ∝ 2^n/2Ω^-1/2exp( n/(2B_⊥^2)). The scaling exponent of  with n differs from that in Grover's algorithm by a factor ∝ B_⊥^-2, which can be made small with large transverse fields n ≫ B_⊥^2 ≫ 1. 
    
    Crucial distinctions between this case and the Hamiltonian in the  analogue version of Grover's algorithm <cit.>  for the case of multiple target states are the
    non-integrability of our model, and the delocalized nature of the
    eigenstates within the energy band W.  Furthermore, analogue 
    Grover's algorithm for multiple targets is  exponentially sensitive in n to the
    weight of the driver Hamiltonian,  and cannot be initialized with a
    computational basis state. 
    
    The model  (<ref>) considered in the paper belongs to the class of n-local infinite range spin glasses similar to   quantum Random Energy Model in transverse field <cit.>. However the key feature of our analysis  –transport via miniband of non-ergodic delocalized states at the tail of the density of states dominated by deep local minima – 
    is ubiquitous to a  broad class  of quantum spin glass models (<ref>), such  as transverse field Sherrington Kirkpatrick, p-spin model <cit.>, K-Satisfiability, etc.
    
    In the above models one can identify    two distinct energy scales. The first scale is the  typical change in classical energy  corresponding to one bit flip: _ flip. The second scale is the typical width of  non-ergodic minibands Γ < Δ_ cl, which decreases   exponentially with n.
    The tunneling transitions between the states inside the  miniband  require  a large number of spin flips, and therefore 
    _ flip≫Γ. Starting from the initial state |z_i⟩ inside the strip of energies Δ_ cl,  the quantum evolution   is confined within the corresponding miniband.  The quantum  PT can be described by an effective down-folded Hamiltonian _ij 
    defined over a subset of computational basis states whose classical energies lie within the energy strip Δ_ cl at the tail of the density  of states. 
    
    We note that once a computational  problem contains a structure, the associated   minibands   can be organized in a   more complex way than in the IB model  considered in our paper. E.g., the population transfer can proceed via the tree of resonances <cit.>. In the structured problems the typical  tunneling matrix elements _ij  can be exponentially greater in n than those in the Grover algorithm  and than  the transition rates in  the classical local search algorithms.
    Extensions of our approach  for the analysis of  the computational complexity of Population Transfer for generic spin glass models  presents a promising direction for  the future research.  
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    Authors  are grateful to Edward Farhi, Lev Ioffe, Vladimir Kravtsov, Christopher Laumann, and  Antonello Scardicchio for the fruitful discussions of this work. K.K. acknowledges support by NASA Academic Mission Services, contract number NNA16BD14C.  This research is based upon work supported in part by the AFRL Information Directorate under grant F4HBKC4162G001 and the Office of the Director of National Intelligence (ODNI) and the Intelligence Advanced Research Projects Activity (IARPA), via IAA 145483. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of ODNI, IARPA, AFRL, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental  purpose notwithstanding any copyright annotation thereon.
    
    .
    
    
    
    IEEEtran
    
    §  MATRIX ELEMENTS OF THE DOWNFOLDED  HAMILTONIAN AND THE NORMALIZATION CONDITION FOR ITS EIGENVECTORS
    
    We introduce eigenstates |x⟩ of the transverse field (driver) Hamiltonian 
    
        H_D=- ∑_j=0^n σ_j^x=∑_x∈H_D^x |x⟩⟨x|,
    
    Here
    
        |x⟩=|x^1⟩⊗…⊗|x^n⟩ ,
    
    where |x^k⟩ is the state of kth qubit  such that σ_x |x^k⟩= (1-2x_k)|x^k⟩ and x-bits take values x^k =0,1. Also 
    
        H_D^x =- (n-2 h_x),   h_x=  ∑_k=1^n x^k,
    
    where h_x is a Hamming weight of the bit-string x and - (n-2 h), h∈(0,n) are eigenvalues of  H_D.
    
    
    We expand  the eigenstates |ψ⟩ of  the system Hamiltonian H (<ref>)  into the basis of the 
     eigenstates |x⟩
        |ψ⟩=∑_x∈Ψ(x)|x⟩ .
    
    We write the  Schrodinger equation H|ψ⟩=E|ψ⟩ in the form
    
        H_D |ψ⟩+∑_j=1^M(z_j)|z_j⟩ψ(z_j) =E |ψ⟩ ,
    
    where ψ(z_j)=⟨z_j|ψ⟩. 
     Then we multiply it  from the  left by ⟨x| and obtain Ψ(x) in terms of ψ(z_j)
        Ψ(x) =∑ _j=1^M (z_j)υ_x,jψ(z_j)/ E- H_D^x .
    
    In Eq. (<ref>) the coefficients  υ_x,j =⟨x|z_j⟩ equal
    
        υ_x,j =2^-n/2(-1)^x· z_j,   x· z_j ≡∑_k=1^nx^k z_j^k,
    
    and  z_j^k=0,1.
    
    We now multiply Eq. (<ref>) from the left by ⟨z_j| where  j ∈(1,M) enumerates marked states and obtain
    
        ∑_x∈ H_D^x Ψ(x)υ_x,j=(E-(z_j))⟨z_j|ψ⟩ .
    
    Plugging here the expression for Ψ(x) (<ref>) the matrix eigenvalue problem (<ref>)  we obtain 
    
        (z_i)ψ(z_i)-∑_j=1^M(z_j) c_i j(E)  ψ(z_j)=E ψ(z_i)
    
    where
    
        c_i j(E)=∑_x∈υ_x,iυ_x,jH_D^x/E- H_D^x,
    
    Because H_D^x depends on a bit-string x only via its Hamming weight ∑_j=1^n x^j one can perform the partial summation in (<ref>) getting
    
        c_i j(E)≡ c(E,|z_i-z_j|),   |z_i-z_j|=∑_k=1^n|z_i^k-z_j^k |,
    
    where the function c(E,d) has he form
    
        c(E,d)=∑_k=0^n-d∑_l=0^dnkn-dl(-1)^l  2^-n/1+E/(n-2 k-2 l).
    
    Above |z_i-z_j| denotes the Hamming distance between the bit-strings z_i and z_j.
    We introduce the rescaling 
    
        ψ(z_i)=A_i/√((z_i)),   i∈[1..M] .
    
    Then Eq.  (<ref>) can be written in the form
    
        ∑_j=1^M_ij(E)A_j=E A_i,
    
    where  _ij is a symmetric M× M matrix
    
        _ij(E)=δ_kj(z_i) +√((z_i)(z_j)) c(E, d_ij),
    
    indices k,j=1:M and δ_kj is Kronecker delta. This is a nonlinear eigenproblem given in the main  text, Eq. (<ref>).
    
    
    
    
    We note that  the projections of the eigenvectors  of  H onto  the  marked state subspace are not, in general, normalized nor they are orthogonal.
    Let us consider the eigenstate |ψ_β⟩ and the corresponding eigenvalue E_β of H. We calculate the corresponding amplitude Ψ_β(x)  using Eq. (<ref>) and plug it into the  normalization condition
    
        ∑_x∈Ψ_β^2(x)=1 ,
    
    obtaining after  partial summation 
     
        ∑_i,j=1^M_i_j  r(E_β,d_ij)ψ_β(z_i)ψ_β(z_j)=1 ,
    
     where the coefficient r(E,d) equals
     
        r(E,d)=2^-n∑_k=0^n-d∑_l=0^d (-1)^kdkn-dl/(n-2(k+l)+E)^2 .
    
     It can be written in the form
     
        r(E,d)=∂/∂ E(c(E,d)-δ_d,0/E) ,
    
     where δ_d,0 is the Kronnecker delta. We use (<ref>) and write
     
        r(E,d_ij)=1/√(_i_j)∂_ij(E)/∂ E .
     
     We now define  the coefficients Q_ij(E) such that
      
        1/ Q_ij(E)=_i_j  r(E,d_ij)=√(_i_j) ∂_ij(E)/∂ E .
    
    Then Eq.(<ref>) takes the form
    
        ∑_i,j1/Q_ij(E)ψ_β(z_i)ψ_β(z_j)=1 .
    
     The above equations   (<ref>) and  (<ref>)  correspond to  Eqs. (<ref>) and (<ref>) of the main text.
     
    
     
    
    
    §  DETAILS OF THE WKB ANALYSIS OF THE COUPLING COEFFICIENTS
    
    
     In the main text we expressed the  coupling coefficient c(E,d)  in terms of the  off-diagonal matrix elements of the resolvent (<ref>) of the transverse field Hamiltonian  H_D between the states that belong to  a maximum total spin subspace S=n/2. The results are given in 
      the expressions (<ref>), 
     (<ref>) from the main text repeated below for convenience
     
        c(E,d)=δ_d,0-E/√(nd) G_n/2 -d,n/2(E)     .
    
    Here the resolvent G_n/2 -d,n/2(E) obeys the inhomogeneous  equation
     
        δ_m,n/2+∑_s=± 1u(m-s/2)  G_m+s,n/2=E G_m,n/2,
    
    
        u(m)=-√(L^2-m^2),   L=n+1/2.
    
    
    We will solve the above equations  for the  case where the  energy E of the resolvent is not far from the center of the Impurity Band 
    
        E=-n+Δ,  Δ=(n^0) .
    
    The WKB solution to Eq. (<ref>)   is sought in the exponential form 
    
        G_m,n/2∝exp(i ∫^m dk  p(k)  ).
    
    It is assumed that ∫^m_0 dk  p(k) = O(n) and |p(k)|= O(n^0) so that G_m,n/2 is varying  steeply with  m changing  by 1.  However  |p^'(m)|= O(1/n) and p(m) is varying very slowly with m  due to the similar property of the coefficients u(m)/L in the Eq. (<ref>). This property is at the root of WKB approximation <cit.>.
    The quantity p corresponds to the  momentum" of the 
     effective mechanical  system with coordinate m, energy E and Hamiltonian function
     u(m)cos p. The function p=p(E,m) is obtained from the equation
     
        u(m)cos p=E.
    
     This equation also defines the  curve on the (m,E) plane with  p=0 shown in Fig. <ref>. Points on that curve are turning points of the classical motion   with  energy E.
     
     For  not too small transverse fields 
        >2L/|E|≃ 1 ,
    
    the Eq. (<ref>)  has two types of WKB solutions that correspond to   real or  imaginary momentum p(m) depending on the value of m relative to the  turning points m=± m_0(E) given below
    [
    The expression (<ref>) for m_0(E) should be evaluated for
    the energy  E=E^0)≃ -n-B_⊥^2 corresponding to the  eigenvalues of impurity band (<ref>),(<ref>)
    In this case the O(n^0) corrections in   the r.h.s. of  (<ref>) vanishes. 
    
        m_0(E^(0))≃n/2√(1-B_⊥^-2)+1-3/4n√(1-B_⊥^-2) .
    
    The terms  O(1/n) above and can be neglected in the exponents of the WKB solutions (<ref>) and (<ref>).]
        m_0=√(L^2-(E/2)^2) .
    
      
    
    
     In the region
     
        n/2+m_0> d> n/2-m_0,
    
     the amplitude    G_n/2-d,n/2 (<ref>) is rapidly oscillating with d and can be written in the form
    
     
     
     
        G_n/2-d,n/2 =-(E)sinϕ(E,d)/ [m_0^2(E)-(n/2-d)^2]^1/4 ,
    
     where
     
        ϕ(E,d)=∫_n/2-d^m_0dk arcsin(√(m_0^2-k^2/L^2-k^2))-π/4 ,
    
    
    
     is a phase of WKB solution and  (E) is the constant of integration that will be discussed  below.
     
     
     On the other hand, in the two regions 
     
        d∈[0,n/2-m_0]∪ [n/2+m_0,n]
    
     the resolvent  G_n/2-d,n/2  is decreasing exponentially with d.
     For example,   in the  left  region 
    
        G_n/2-d,n/2= (E)/2e^ | Im ϕ(E,d)|/[(n/2-d)^2-m_0^2(E)]^1/4 .
    
    
    
    We omit here for brevity the expression in the right  region (<ref>).
    
    
    
     §.§ Determination of the integration constant in WKB solution
    
    Within the  WKB approach the  integration constant   (E)  can be obtained by matching the exponential asymptotic  (<ref>) with the 
    solution obtained near the boundary of the interval d=0. However as discussed in Sec. <ref> of the main text,  for the  relevant range of the model parameters
    the properties of the typical sample in  the ensemble of the IB Hamiltonians  depend only on G_n/2-d,n/2 in  the  region of its oscillatory behavior (<ref>) away from the boundaries of the interval d=0, n. To avoid the analysis in the region of no consequence for us  we determine (E) by equating  
    the above WKB asymptotic for  G_n/2-d,n/2 at the center of the interval d=n/2 with  expression for G_0,n/2  at that point obtained  in a different way.
    
    Using  Eq. (<ref>) we  write c(E,n/2) in the integral form
    
        c(E,n/2)= i E/2^n  ∫_0^∞dτ (1-e^4iτ)^n/2 e^i(E/-n+i o)τ
    
    (o→ +0). The integral can be expressed in terms of the Gamma function Γ(x). In the region of not too small  transverse fields (<ref>)
    it has the form
    
        c(E,n/2)=2^1-n π a (a^2-1)^-1Γ( n/2 )  /sin( π(a-1)n/4a) Γ((a+1)n/4a)Γ((a-1)n/4a).
    
    where 
    
        a=-n /E >1 .
    
    Using Sterling  formulae for Gamma function we obtain in the limit n≫ 1, a=(n^0)
        c(E,n/2) =√(nπ)/2 asin( π(a-1)n/4a)2^-n/2 e^-n θ(a),
        
           θ(a) =2arctanh(1/a)+aln1-a^-2/4a .
    
    For  large transverse fields a≫ 1 and we have  θ≃ a^2/4.
    
    Using Eq. (<ref>) we obtain the asymptotic of the Green function at the zone center
    
        G_0,n/2(E)=π/8 n^3^1/4exp(-n θ(a))/√(-1) sinϕ(n/2,E) .
    
    Here we used the equality for the phase WKB ϕ(E,n/2) (<ref>) at the zone center
    
        ϕ(E,n/2)=π(a-1)n/4a .
    
    On the other  hand, from the WKB expression (<ref>)
    we get
     
        G_0,n/2 =-(E) 2/n^1/2 sinϕ(E,n/2)/ (1-B_⊥^-2)^1/4 .
    
    By comparing the Eqs. (<ref>) and (<ref>) we finally obtain the constant of integration (E)
        (E)=-π^1/4/(32n (-1)^1/4exp(-nθ(a))/(sinϕ(E,n/2))^2.
    
    
    
    
    One can use   (<ref>) in (<ref>) and (<ref>)  to obtain the  expression for c(E,d)  in the region (<ref>). Before providing  the result 
    we observe that for energies E  not too far from the Impurity Band center (cf. Eq. (<ref>)) the expression for nθ(a) can be expanded in powers of 1/n
        nθ(a)≃ nθ()-E+n/2+(n^-1) .
    
    where E+n≡Δ=(n^0).
    
    Finally the expression for the coupling coefficient has the form
    
        c(E,d) =   √(A(E,d/n))  n^1/4 e^-nθ()/√(nd)
        ×   √(2) sinϕ(E,d) ,
    
    where  the WKB phase ϕ(E,d) is given in (<ref>) and the coefficient A(E,ρ) equals
    
        A(E,ρ)=√(π/32 )e^E+n//(B_⊥^2-1)υ(ρ)sin^4(ϕ(E,n/2)) ,
    
        υ(ρ)=(1-(1-2ρ)^2/1-B_⊥^-2)^1/2 ,
    
    It is related to  A(ρ) in the Eq. (<ref>) of the main text as follows: A(ρ)= A(E^(0),ρ).
    The phase  ϕ( E,n/2) in (<ref>) has an explicit form
    
        ϕ( E,n/2)=π/4(n(1-B_⊥^-1)+n+E/) .
    
    
     §.§ Limit of large transverse fields ≫1
    
    In the limit of large transverse fields the tuning point m_0 (<ref>) is very close to the boundary of the interval m=L so that one has a small parameter 
    
        √(L-m_0/L)=1/√(8)|E|/L≪ 1
    
    In this case the expression for the WKB phase takes a simple  form
    
        ϕ(E,d)=π d/2-π n/4χ(E,d/n)/ ,
    
    
        χ(E,ρ)= 1-Δ/n1-2/πtan^-11-2ρ/√(1-(1-2ρ)^2)
    
    where Δ=E+n=(n^0) and values of d are not too close to the interval boundaries
    
        n-d,  d≫ L-m_0∼n/ .
    
    We note that for large transverse fields ≫1 the phase is a sum of the two terms. First term changes rapidly with d with the slope π/2 and second term changes very little (by an amount (n^-1)) when d is changed  by 1. 
    
    
    
    We note that unlike the study of the WKB eigenfunctions  where one has to  select the WKB solution that  decays into the classically forbidden region (<ref>), 
    the Green function  G_n/2-d,n/2(E) corresponds to the  solution that  increases exponentially with  m=n/2-d>m_0.  
    Using the oscillating (<ref>) and exponentially growing  (<ref>)  WKB solutions   one can obtain the  coefficient c(E,d) from the relation (<ref>).
    
    
    
    This  will provide an asymptotic  WKB form of c(E,d)   almost everywhere on the interval d∈ [0,n] except  for the small vicinities  of the turning  points, | n/2-m_0(E)-d |= O(n^0) and end points, n-d, d= O(n^0).
    In Fig. <ref> we plot the comparison between the coefficients c(E,d) computed based on exact expression (<ref>) and the results of asymptotic  WKB analysis using Eqs.  (<ref>),(<ref>).  
    
    
    
    
    
    § LINEARIZATION OF THE DOWN-FOLDED HAMILTONIAN NEAR THE CENTER OF THE IMPURITY BAND
    
    
    We divide   the  Hamiltonian H(E) for a given E  on two parts, accordingly
    
        _ij(E)=^(0)_ij(E)+_ij^(1)(E) ,
    
    where we defined
    
        ^(0)_ij(E)=n(c(E,0)-1) δ_ij,
    
    
        _ij^(1)(E)=δ_ij(1-c(E,0))ϵ_i+n c(E,d_ij)(1-δ_ij).
    
    We write similar expansions for energies and amplitudes
    
        E≈ E^(0)+E^(1),  ψ(z_j)≈ψ^(0)(z_j)+ψ^(1)(z_j),
    
    and get
    
        (E)  ≈^(0)(E^(0))+^(0)(E^(0))/∂ EE^(1)+^(1)(E^(0)),
    
    where the parts of the Hamiltonian ^(0,1) are given above.
    We plug the above expansions  into the system of equations  (<ref>)  ∑_j=1^M_ij(E)_j=E _i, and  use   (<ref>) to express   _j^(0)=n^1/2ψ^(0)(z_j). 
    Equating terms of the same order in ϵ_j and c(E,d_ij), i≠ j, we obtain the equation for eigenstates and eigenvalues in  zeroth order
    
        n[c(E^(0),0)-1]ψ^(0)(z_j)=E^(0)ψ^(0)(z_j) ,
    j∈[1..M], and in the first order
     
        a ϵ_i ψ^(0)(z_j)+b∑_j≠ i=1^Mn  c(E^(0),d_i j)ψ^(0)(z_j)
        
        =E^(1)ψ^(0)(z_j)
    
    Above index j enumerates marked states. Also the coefficients a,b equal
    
        a=b (1-c(E^(0),0),   b^-1=1-n∂ c(E^(0),0)/∂ E .
    
    
    
    
    Similarly to the above we find from Eqs. (<ref>),(<ref>) the zeroth-order approximation to the the total probabilistic weight  of an eigenfunctions |ψ⟩  over the marked state  subspace Q_jk^(0)=δ_jkQ where
    
        Q_jk^(0)=δ_jkQ,  1/Q =n^2 ∂/∂ E(c(E,0)-1/E)_E=E^(0).
    
    
      §.§.§  Zeroth-order of the perturbation theory
    
    
    Eq. (<ref>) admits the solution corresponding to the M-fold  degenerate energy level that originates from the band of the  marked  states, E^(0)→ -n
    in the limit of → 0. The  corresponding M eigenstates ψ_β(z_j) (β∈[1..M]) have support  over the  part of  computational basis   corresponding  to  
    marked  states:   ψ_z_j^β≠ 0,  j∈(1,M).  
    Using c(E,0) from (<ref>) the explicit form of the equation (<ref>) for eigenvalue in zeroth order is given in the main text, Eqs. (<ref>),(<ref>) which we repeat here for convenience. 
    
        E^(0) =-n-Δ_0 ,
    
    
        Δ_0=n 2^-n∑_d=0^nnd(n-2d)/n+Δ_0-(n-2d)  .
     
     Here    Δ_0 is  the root of the  above  transcendental equation that satisfies the condition  lim_→ 0Δ_0=0. In general,  
      the sum (<ref>) is dominated by the region of   values of d such that |d-n/2|= O(n^1/2). We    obtain   Δ_0  in a form of a series expansion in powers of n^-1
        Δ_0≃ -B_⊥^2-B_⊥^4/n+ O(n^-2),
      
     
    
    Similarly, using  c(E,0) from (<ref>) in the equation (<ref>) for the zeroth-order total weight over the marked state subspace we obtain
    
        ∑_k=1^M|ψ_z_k^(0)|^2=Q,
    
    
        1/Q=1/2^n∑_d=0^nnd1/((n-2d)-n-Δ_0)^2.
     
    Using (<ref>) and employing   similar approximations to that from the above we get an asymptotical expression in large n limit
    
        Q≃ 1-B_⊥^2/n-3B_⊥^4/n^2+ O(n^-3).
    
    We recall that in our study n is asymptotically large and  we  always  assume  that the transverse field = O(n^0) (but can be parametrically large, ≫ 1).
    
    
    The denominator in Eqs. (<ref>),(<ref>) corresponding to d=m will become zero  at  resonant" transverse field value =B_⊥  m   which is a root of the equation (<ref>) in the main text.
    
    
    
      In the range of  under consideration n/2-m≫ n^1/2. 
    
    Near the mth resonance the term with d=m in the sum (<ref>) becomes anomalously large due to a small denominator despite the factor p_m being very small.
    We keep this term  (<ref>) along with the   terms corresponding to |n/2-d|∼ n^1/2 and obtain 
    
        Δ_0 ≃δ_B/2±√(δ^2_B/4 +n^2 p_m) ,
    
    where we  introduced  rescaled transverse field difference from its value at resonance
    
        δ_B=n-B_⊥  m/B_⊥  m^(0) ,
    
    where B_⊥  m^(0)=n/(n-2m).
    
    
    Clearly, in the resonance region  δ_B∼  n  p_m^1/2 and |-B_⊥  m|∼Δ B_⊥  m where Δ B_⊥  m∼  2^-n/2nm B_⊥  m^(0).
    There   the weight factor Q is decreasing dramatically (cf.   Fig. <ref>) and   the above perturbation theory breaks down.  The width of the resonant regions  Δ B_⊥  m   (<ref>)   remains exponentially small in n  for   n/2-m≫ n^1/2.
    
     In this study we will only focus on the off-resonance case, assuming the condition
    
        Δ B_⊥  m≪  |B_⊥  m- |∼ | B_⊥  m+1- | = O() .
    
    
      §.§.§  First order of the perturbation theory
    
      *  In the fist order of perturbation theory we obtained effective Hamiltonian linearized in the vicinity of shifted band.
    
      *  give the expression for this Hamiltonian: only introduce terms that depend on n and explain. Temove θ symbol given assymtotics.
    Give cartoon of the of the levels  of the driver and impurity band.
    
    The first order equation (<ref>) determines the correct zeroth order eigenstates {ψ_β(z_j)}_β=1^M and removes the degeneracy of the energy levels. To evaluate the coefficients a,b in (<ref>) we calculate c(E,0) away from resonance using the same approach as that in the evaluation of the sum
    in (<ref>)
    
        c(E,0)≃ -n B_⊥^2/E^2+ O(n^2 B_⊥^4/E^4 ) .
    
    The coefficients    a,b≃ 1+ O(B_⊥^2/n) and in what following will be  replaced  by unity. Then Eq.  (<ref>) corresponds to the effective Hamiltonian  with the  matrix elements, _ii=_i and _ij≠ i=n c(E^(0),d_ij) where coupling coefficients c are  given in (<ref>). 
    
    
    
    Using Eqs. (<ref>),(<ref>) for zeroth-order energy E^(0),
    the matrix _ij  can be written in the form (<ref>). 
    
    
    
    
    § STATISTICAL INDEPENDENCE OF MATRIX ELEMENTS
    
    In this paper the IB  Hamiltonian _ij is determined by the symmetric matrix of Hamming distances d_ij between the bit-strings corresponding to the marked  states sampled without replacement   from the set of all possible  2^n bit-strings.
    
     Instead of this  ensemble one can consider   a different one,  where each of the M bit-strings   is sampled  with replacement from the full set . In this ensemble Hamming distances  d_ij for distinct pairs i,j  are statistically independent allowing for much simpler statistical averaging. Indeed, for a given row i of the  matrix  d_ij  the  joint probability distribution of the  two distinct off-diagonal matrix elements can be estimated as,
    
        p_d_i j_1,d_i j_2-p_d_i j_1 p_d_i j_2∝1/2^nΔ(d_i j_1-d_i j_2) p_d_i j_1.
    
    Here Δ(d) denotes the Kronecker delta,   j_1≠ j_2≠ i and p_d as before corresponds to the modified binomial distribution,
      
        p_d =1/Z 2^-nnd,   Z=∑_d=1^n2^-nnd,
    
     (also   ∑_d_1,d_2=1^np_d_1,d_2=1).
     One can see that  the statistical correlation between a pair of Hamming distances d_i j_1,d_i j_2  is exponentially small (in n) and can be neglected.
      
     Such an ensemble allows for multiple copies of the same bit-string to be sampled. However this effect is not statistically significant for  modest values of M
        1≪ M  ≪ 2^n/2 .
    
     This can be seen by comparing the number of ways to perform unordered sampling of M elements from the group of 2^n elements with and without replacement. Using Stirling's formula we write the former number as,
     
        2^n+M-1M≃2^nM exp(M^2/2^n)(1+ε),
    
    where the latter number is given by 2^nM with ε∼ M 2^-3n/4≪ 1.  It is clear that when condition (<ref>) is satisfied the two ensembles are statistically equivalent because repetitions  can be neglected.
    
    
     
    
    §  BOUND ON THE LARGEST EIGENVALUE OF _IJ FROM GERSCHGORIN  CIRCLE THEOREM
    
     
    One can use the above estimates of the typical largest matrix elements of the matrix _ij    to consider  the  bounds on its   eigenvalues  given by the  Gerschgorin  circle theorem <cit.>.
      For the  case of real  eigenvalues the theorem states that every eigenvalue   lies within at least one of the intervals [_ii- R_i,_ii +R_i ] where i∈[1..M] and R_i=∑_j≠ i|_ij| is a sum of absolute values of the off-diagonal elements  in the ith row.  
      For a  randomly chosen row the value of  R_i can be estimated as follows 
    
        R_i ≃ M∑_d=1^n  p_d |V(d)| ,
     
    where p_d is defined  in (<ref>). From Eq. (<ref>) one can see that the above sum is dominated by the terms satisfying |n/2-d|≪ n. Using Stirling's approximation  we get 
    
    R_i ∼ M 2^-n/2 e^-nθ.  
    
    For typical diagonal matrix elements |_ii|=|ϵ_i| W. 
     Therefore from the Gerschgorin  theorem we conclude the eigenvalues E_β^(1) of   satisfy the following bound 
     
        | E_β^(1)| ≤max{W,  M 2^-n/2 e^-nθ} .
    
    One can see that Gerschgorin bound in our case precisely corresponds to the typical maximum  element in the matrix _ij.
    
    
    
    
    
    
    §  MEAN VALUE AND STANDARD DEVIATION OF THE OFF-DIAGONAL MATRIX ELEMENTS _IJ
    
     
     The mean value of the off-diagonal matrix element 
    
        ⟨_ij⟩=n∑_d=0^np_d c(E,d)≃n/2^n/ -  1.
    
    is much smaller than its  standard deviation
    
        ⟨ ( _ij-⟨_ij⟩)^2⟩^1/2≃√(n/2^n).
    
    This is related to the symmetry p_d=p_n-d and a rapid oscillation of   c(E^(0),d) with d (cf. (<ref>), (<ref>) and Fig. <ref>). 
    
    
    
    We note from (<ref>),(<ref>) that the standard deviation is exponentially larger than the typical value 
    
        ⟨ (_ij-⟨_ij⟩)^2⟩^1/2∼  V_ typ e^n θ.
    
     This can be understood by looking at the  values of d_ij that dominate the  variance of _ij.
     We write
    
        ⟨ (_ij)^2⟩ = n^2∑_d=0^n c^2(E^(0),d) p_d.
    
    It follows from the Eqs.  (<ref>) and   (<ref>) that   for   d∈(n/2-m_0,n/2+m_0) the coefficient c^2(E,d)∝ 1/nd decreases exponentially with d, while
    the distribution p_d ∝nd increases exponentially with d. The binomial factors cancel out and  the expression under the summation in (<ref>)  
    contains very slowly-varying with d (non-oscillatory) part. However for  d∈(0,n/2-m_0)    the coefficient c(E,d) grows exponentially  faster than 1/nd with decreasing d (see Eqs. (<ref>), (<ref>)).   Therefore the  variance (<ref>) is dominated by  non-extensive values of d= O(n^0)
    that are  much smaller than the smallest Hamming distance d_ min= O(n)  (<ref>) in a  randomly chosen  row of d_ij. Therefore the variance of _ij is not a good statistical characteristic of the PDF  of _ij. It  is dominated by the extremely rare atypical instances of the ensemble.
    
    
    
    
    
    
    
    
    
    § PDF OF THE SQUARED OFF-DIAGONAL MATRIX ELEMENTS OF IMPURITY BAND HAMILTONIAN
    
    
    
    In this section we provide the details of the derivation of the PDF for the non-oscillatory parts of the (squared) off-diagonal matrix elements _ij^2 of the IB Hamiltonian. As discussed in the main text,  in the  asymptotical limit of large n≫ 1 one can make an approximation that  n  is a continuous variable and  we replace the summation over d in (<ref>) by an integral and Kronecker delta 
     δ(x) by Dirac delta. This results in the  Eq. (<ref>) displayed below for convenience
     
        P(_ij^2)= ∫_0^n p_x   δ(V^2(x)-_ij^2)dx .
    
    It was discussed in the main text (see also below)  that the condition for this validity of this approximation is
    
        1/nlog_2 M≪ 1 .
    
    It corresponds to the number of marked states M that is not very large. For example, it can still scale exponentially with n so that M=2^μ n, μ= O(n^0),
    but the coefficient μ in the exponent needs to be small μ≪1. 
    
    
     The  expression (<ref>) is obtained  using the analytical continuation p_x of the binomial distribution  p_d (<ref>)  from the integer domain d∈[0,n]  onto  the interval of a real axis x∈[0,n] in terms of the Beta function B(x,y)
        p_x=2^-nnx=2^-n/(n+1)B(x+1,n+1-x) ,
    
     and the resulting identity
    
        ∫_0^ndx  p_x=1.
    
    
    In what following we will study the   rescaled quantities 
    
        w_ij≡_ij^2/V_ typ^2=(2/π n)^1/21/p_d_ij ,
    
    where  i≠ j,  V_ typ  is  given in (<ref>) and p_d=2^-nnd. Using Stirling's approximation in binomial coefficient
    
        p_x≡ p_B(x/n),   p_B(ρ)=e^-n (ρ)/√(2π nρ(1-ρ)) .
    
    
        (ρ)=ρlogρ+(1-ρ)log(1-ρ)+log 2 .
    
    we get from Eq. (<ref>) for _ij=V(d_ij)
        (ρ) ≡V^2(n ρ)/V_ typ^2≃√(4ρ(1-ρ))/υ(ρ) e^n (ρ) ,
    
    where   υ(ρ) is given in (<ref>).
     Eq. (<ref>) takes the form
    
        w_ij=(d_ij/n) .
    
    Then the expression for the PDF for w_ij
        g(w_ij)=V_ typ^2 P(V_ typ^2 w_ij) ,
    
    can be written in the form (cf. (<ref>))
    
        g(w)=2n∫_0^1/2p_B(ρ)δ(w-(ρ))dρ .
    
    We note that the domain of g(w) is bounded from below by w=1 and from above by w= O(2^n).  Taking the integral in (<ref>)  we get
    
        g(w)=2np_B(ρ_w)/ | d(ρ)/dρ |_ρ=ρ_w ,
    
    where   the rescaled Hamming distance ρ_w  is a root of the transcendental equation
    
        (ρ_w)=w .
    
    In the leading order in n≫ 1 this equation gives
    
        (ρ_w)=1/nlog w ,
    
    where (ρ) is given in (<ref>).
    Also  using Eqs. (<ref>),(<ref>) in (<ref>) we  get 
     
        g(w)=1/w^2√(πℓ(w)) ,
    
    where 
    
        ℓ(w)=n/8υ^2(ρ_w) |log(ρ^-1_w-1)|^2
    
    Here   the dependence  of  ℓ(w) on   w is shown in Fig. <ref>. In the entire range the dependence is logarithmically  slow. 
    
    We note that the (<ref>) is a valid approximation to (<ref>) for ρ-ρ_0≫ 1/n where ρ_0 is a zero of υ(ρ)
        υ(ρ_0)=0,  ρ_0=1/1-√(1-B_⊥^-2)
    
    It corresponds to Hamming distance n ρ_0 = n/2-m_0  (<ref>) which lies at the boundary of the interval  (<ref>) where the WKB solution (<ref>),(<ref>) applies (see the discussion in the Sec. <ref>). It is assumed that n ρ_0 is smaller than the typical smallest Hamming distance  in a randomly selected row
    
        - n ρ_0≫ 1.
    
    Using the asymptotic expression (<ref>)  for the binomial distribution in Eqs. (<ref>)  we get the equation for  in the form
    
        (/n)=1/nlog M
    
    The function (ρ) is decreasing with ρ for ρ(0,1/2). Therefore Eq. (<ref>) leads to  the condition (ρ_0)-(/n)≫ 1/n, or
    
        (ρ_0)-1/nlog M≫1/n .
    
    Using explicit forms of (ρ) and ρ_0 we get in the limit of ≫ 1
        log 2-1/nlog M     >2log+2log 2+1/4B_⊥^2+ε
        
        0<ε   = O(B_⊥^-4) .
    
    This is the condition for (<ref>). Clearly it corresponds to  a much weaker constraint on the values of M  than the condition 1/nlog M<1/2log 2
    provided by the requirement of a statistical independence of matrix elements of _ij (cf. (<ref>)).
    
    
    
    
     
    
      §.§.§  Case of 1/nlog_2 M ≪ 1
    
    
     The rescaled Hamming distance ρ_w  depends on w via the logarithmic factor α=1/nlog_2w. This dependence is shown in the 
      inset to the Fig. <ref>.  In this section we consider  
      
        α=1/nlog_2 w≪ 1 .
    
      Then we get
      
        p_B(ρ)≃(2/π n)^1/2 e^-n (ρ) ,
    
    
        (ρ)≃ 2(1/2-ρ)^2 .
    
      Then using (<ref>) we get
      
        ρ_w ≃1/2-(α/2)^1/2 .
    
    
        ℓ(w)≃log w ,
    
      and finally,
      
        g(w) ≃ g_∞(w)=1/w^2 √(πlog w) .
    
    The subscript here indicates that, unlike g(w), the PDF g_∞(w) has the upper boundary of its domain  equal to infinity.
    It is of interest to calculate  for a given w the magnitude of the relative changes of   V^2(d) and of the binomial coefficient p_d when Hamming distance is chaining by 1 (and ρ_w is changing by 1/n). We define as in Eqs. (<ref>),(<ref>) w=V^2(d)/V_ typ^2 and obtain
    
        V^2(d+1)-V^2(d)/V^2(d)≃p_d+1-p_d/p_d
           ≃ 4(1/2-ρ_w )=√(8α)≪ 1.
    
    Here we used Eqs. (<ref>) and (<ref>).  The above inequality justifies using the continues approximation (<ref>) in (<ref>).
    
    
    
     In a randomly chosen row of w_ij  the PDF that  the  largest element equals w is
     
        PDF(max_mw_m=w)≃M e^-M wlog w/wlog w,   M≫ 1 .
    
    Typical largest  element in a row max_1<j<iw_ij∼ M in agrement with the results obtained earlier,  cf. Eqs. (<ref>), (<ref>) and (<ref>). 
    Therefore in order to ensure that  α≪ 1 for all matrix elements in a typical row of w_ij we require that log_2 M ≪ n
        1≤ w≲  M,  1/nlog_2 M≪ 1 .
    
    
    
    The typical value of the smallest element in a randomly selected row of the rescaled matrix of Hamming distances d_ij/n equals
    
        ρ_ min=/n=1/2-√(log_2 M/2n) .
    
    We note that in the case we consider 
    
        n ≫ n/2- =  O(n) ,
    
    minimum value  is close to n/2 but is still separated by extensive distance from it.
    
    In this paper we use the  expression for  the  matrix elements of the IB Hamiltonian  _ij (<ref>) that only applies in the region |n/2-d_ij|<m_0
     where  m_0  is given in (<ref>).  The elements in a typical row of the matrix d_ij belong to this region if the condition |n/2-|<m_0 is fulfilled.
     Using (<ref>) we can re-write this as an inequality for M
        M<2^n/2(1-B_⊥^-2) .
    
     This inequality is satisfied under the condition (<ref>).
     
    
    
    
    
    
    
    
    §  CHARACTERISTIC FUNCTION OF THE PDF OF  THE SQUARED OFF-DIAGONAL MATRIX ELEMENTS OF IMPURITY BAND HAMILTONIAN
    
    
    Here we compute the characteristic  function of the PDF  g_∞(w) (<ref>) (also given in Eq. (<ref>) of the main text). It is defined as follows
    
        ϕ_∞(u)=∫_1^∞dw  g_∞(w)(e^i u w-1) .
    
    We will be interested in the asymptotic limi of the above expression at small |u|≪1. It is convenient  to calculate separately real  and imaginary parts of ϕ_∞(u). 
    
    
    For real part we have
    
        -√(π)/2  Re[ϕ_∞(u)]=∫_1^∞1/x^2√(log x) sin^2(u x/2) .
    
    Because ϕ_∞(-u)= ϕ_∞^*(u) we can assume that u>0 and break the  interval of integration  above on two parts
    
        x∈ [1,X/u] ∪ [X/u,∞),   u≪ X≪ 1 .
    
    We write
    
        -√(π)/2  Re[ϕ_∞(u)]=R_1(u)+R_2(u) .
    
    Here
    
        R_1(u)=∫_1^X/u1/x^2√(log x)sin^2(u x/2) ,
    
        R_2(u)=∫_X/u^∞1/x^2√(log x)sin^2(u x/2) .
    
    Using (<ref>) asymptotic expansion of R_1(u) has the form
    
        R_1(u)≃u  X/4(log(1/u))^1/2+u  Xlog(1/X)/8 (log(1/u))^3/2 +…
    
    Also after some tedious calculations we obtain
    
        R_2(u)    ≃u/(log(1/u))^1/2(π/4-X/4
        
        +   u/2(log(1/u))^3/2π (γ_ Euler-1)/4 .
     
    where
    
        γ_ Euler≃ 0.577
     
    is the  Euler constant.
    
    
    Similarly to the above we also  break the interval of integration in the imaginary part of ϕ_∞(u) on two parts given in (<ref>)
    
        [ϕ_∞(u)]=I_1(u)+I_2(u) .
    
    where
    
        I_1(u)=∫_1^X/usin u x/x^2√(πlog x) ,
    
        I_2(u)=∫^∞_X/usin u x/x^2√(πlog x) .
    
    Expanding the integrand (<ref>) in u and using condition  (<ref>)  we get
    
        I_1(u)≃2u√(log1/|u|)/√(π)-u log1/X/√(πlog1/|u|)+u log^2X/log^3/2|u| .
    
    Performing similar asymptotic expansion in I_2(u) we obtain
    
        I_2(u)≃u(1-γ_ Euler-log X)/√(πlog1/|u|)+u log^2X/log^3/2u .
    
    
    
    Finally, we combine together Eqs. (<ref>),(<ref>) into Eq. (<ref>) to obtain first two terms in the asymptotic expansion  of Re[ϕ_∞(u)]  in powers of 1/log u≪ 1
        Re[ϕ_∞(u)]=-|u| √(π)/2√(log |u|^-1)1-1- γ_ Euler/2log |u|^-1 .
    
    We also combine together the above expressions for  I_1 and I_2 to obtain a similar asymptotic expansion of  Im[ϕ_∞(u)]
        [ϕ_∞(u)]≃2u√(log1/|u|)/√(π)+u(1-  γ_ Euler)/√(πlog1/|u|) .
    
    Note that in both cases the terms involving X cancels out confirming the validity of the matching procedures.
    
    
    
    §  GENERALIZED CENTRAL LIMIT THEOREM FOR THE SUM OF  M RANDOM VARIABLES W_M THAT OBEY   THE DISTRIBUTION  G_∞(W)
    
    
    In this section we will study the asymptotic PDF  for the sum the independent identically  distributed random variables  in  Eq. (<ref>) sampled from the probability distribution (<ref>).  We note that the variance of the random variables does not exist.
    The  PDFs with polynomial tails at infinity are known as Pareto (heavy-tailed) distributions.
    According to the Generalized Central Limit Theorem (GCLT), the  PDF  of the  sum of M Pareto variables  for M→∞ approaches its asymptotic form   given by the stable law  <cit.>.This general property coincides  with the  usual  Central Limit Theorem for the case when random variables in a  sum have finite variances. In this case the limiting  PDF has Gaussian form. 
    
    We note that  the PDF given by Eq. (<ref>)  is not strictly polynomial at w→∞  because of the additional logarithmic factor. We will  derive the asymptotic form of the sum (<ref>) of random variables (<ref>) explicitly and compare with the standard GCLT result  without the logarithmic factor.
    
    
    We are interested in the PDF of the random variable s_M such that (cf. (<ref>),(<ref>))
    
        s_M =1/M∑_i=1^Mw_i,   g_∞(w)=  1/w^2 √(πlog w) .
    
    Here  w_i are i.i.d random variables sampled from g_∞(w) and we are interested in the asymptotic limit M≫ 1.
     
    Using the convolution property of a sum of statistically independent random variables we get for the PDF   of s_M
        (s_M)=1/2π∫_-∞^∞dq   [φ_∞(q/M)]^M  e^-i  q s_M ,
    
    where 
    
        φ_∞(u)=1+ϕ_∞(u) ,
    
    and ϕ_∞(u) is given in (<ref>).
    The limit M≫ 1 corresponds to |u|≪ 1. We note that
    
        lim_u→ 0ϕ_∞(u)=0
    
    Taking into account that  ϕ(u) is small in the above limit we write
    
        (s_M)≃1/2π∫_-∞^∞dq exp[-i q s_M+Mϕ_∞(q/M)] .
    
    
    
    Quantity  M ϕ_∞(q/M) can be expanded in inverse powers of log M≫1  using asymptotic form of the characteristic function at small argument given in 
    Eqs. (<ref>),(<ref>). First few terms of expansion have the form
    
        M  ϕ_∞(q/M)   ≃ -π |q|/2√(log M)+√(π) |q|  1-γ_ Euler-log|q|/4(log M)^3/2 ,
        
        
        M  ϕ_∞(q/M)    ≃  2q  (log M/π)^1/2+ q 1-γ_ Euler-log |q|/(πlog M)^1/2
           +qlog|q|(1-γ_ Euler ) /2√(π)(log M)^3/2 .
    
    where γ_ Euler is the Euler constant.
    
    
    It is clear from comparing  individual terms in  Eq. (<ref>) with the exponential in the integrand in Eq.(<ref>) that q=(√(log M)). 
    Therefore we can drop in Eqs. (<ref>) terms (1/(log M)^3/2).
    We make the change of variables in the integral in (<ref>)
    
        q=2 √(log M/π) t ,
    
    and obtain
    
    
        (s_M)=1/σ_ML_1^1,1(s_M-b_M/σ_M) ,
    
        L_1^1,1(x)≡1/2π∫_-∞^∞dt e^-i t x-|t|-2 i/πtlog|t| .
    
    Function L_1^1,1(x) above  is  a so-called Levy alpha-stable distribution  <cit.>  shown  in Fig. <ref>. The distribution  is defined by its characteristic function. Parameters b_M and  σ_M in (<ref>) are typical values that  characterize the shift of the maximum of the (s_M) from the origin and its  overall scale, respectively. They are given in 
     Eqs. (<ref>) and (<ref>) of the main text and we also provide them for convenience  below
    
        σ_M=π/21/(πlog M)^1/2 ,
    
    
        b_M≃σ^-1_M-2/πσ_Mlog(σ^-1_M) +2/π (1-γ_ Euler) σ_M ,
    
    where γ_ Euler  is the  Euler constant.
    
    
    It is instructive to compare the above expressions   with the result for the sum 
     of random  variables  that obey a standard Pareto distribution 
     (i.e.,  without the logarithmic factor present in g_∞(w))
    
        s_M^0=1/M∑_i=1^Mw_i,   w_i ∼ g_0(w)=w^-2 .
    
    The PDF of s_M^0 has the same form as the PDF of s_M given in (<ref>)  but the expressions for 
    the shift b_M^0 and  the overall scale σ_M^0 are   different
    
        σ_M^0=π/2,   b_M^0=log M+1-γ_ Euler+log(π/2)
    
    One can see that
    
        σ_M^0/σ_M∼b_M^0/b_M∼  (log M)^1/2≫ 1 .
    
    The rescaling  factor (log M)^1/2 between the PDFs of s_M and s_M^0  can be explained  by a similar logarithmic factor in the ratio g_0(w)/g_∞(w)∼ (log w)^1/2, taking into account the fact that typical of w∼ M.
    
    
    
    
    
    §  JUSTIFICATION OF REPLACING SUM WITH INTEGRAL IN EQ.( <REF>).
    
    
    We note that the number of marked states  in a miniband (<ref>) on a Hamming distance d from a given marked  state |z_j⟩ decreases  rapidly when d.   There is a  typical minimum Hamming distance  d ≃ such that 
    
        d_ min^ res= argmin(Ω_d)=(1) .
    
    There will be  no states in the  miniband located at the Hamming distances
    d from the state |z_j⟩ that lie inside the intervals
     d∈[1,)∪ (n-,n]. For those values of d    we have Γ_j^(d)=0.
    Using  (<ref>) we get
    
        ≃n/2-√(n/2log2A Ω/π n^2) .
    
    where A=A(E^(0),1/2) (<ref>). 
    
    On the one hand we assume throughout the paper that  the number of marked states in a miniband Ω≫1 is sufficiently large so that
    the number n-2 of dominant  terms in the sum    (<ref>) is much bigger than  1.
    For example, using the  scaling ansatz (<ref>) we have Ω∼ M^2-γ (<ref>). Then  assuming that γ<2 and 1>1/nlog_2 M= O(n^0) we can see that 
     the  second term in the r.h.s of  (<ref>) is of the order of n and  therefore the number n-2=(n).
    
    
    On the other hand we note that the number Ω_d (<ref>) of marked states in a miniband  on a Hamming distance d from a given marked  state |z_j⟩  is large (M_j^(d) > ≫1)  for almost all d, aside from  (n^0) values of d near the boundaries of the interval d∈[, n-].  
    
    We recall that  all terms  in a sum  (<ref>) are nearly equal to each other and therefore the relative  contributions to Γ_j from the boundary terms is (1/n) and can be neglected in a leading order estimates of the typical quantities. For d  away from the interval boundaries  the function δ_η(_j-_m) in Eq. (<ref>) changes little between the adjacent values of _m (by an amount ∼ 1/≫ 1). This provides the justification for us to replace  the sum over m in (<ref>) by an integral.
    
    
    
    
    §  PDF  OF  THE  RANDOM VARIABLE H=Η/(Z-)^2+Η^2
    
    
    
    
    Consider the PDF  p_η(h;z)  introduced in the Eq. (<ref>)
    
        p_η(h;z)=∫_-∞^∞1/W p_A(/W) δ[ h-δ(z-,η) ] d ,
    
    Here the function of two arguments δ(x,y) is defined in (<ref>) and  δ[x] is Dirac delta-function denoted here with bold font to distinguish from the above function. We also used the relation (<ref>) for the PDF of marked state energies.
    Solving equation
    
        h=η/(z-)^2+η^2 ,
    
    for  we get
    
        _±=z±√(η(h^-1-η))
    
    From here and from (<ref>) we get
    
        p_η(h; z)=√(η)/2 h^3/2√(1-η h)(φ_+(h;z)+φ_-(h;z))
    
    
        φ_±(h;z)=W^-1p_A(z±√(η(1/h-η)))
    
    For |z|≪ W we get p_η(h;η)≃ p_η(h;0)
        p_η(h;0)=√(η)/h^3/2√(1-η h) p_A(√(η(1/h-η))) .
    
    
      
    Uniform Distribution
    For the case of uniform distribution
    
        p_A()=1/Wθ(W/2-) ,
    
    where θ(x) is Heaviside theta-function we have
    
        p_η(h;0)=1/h^3/2√(η^-1-h) .
    
    Domain of values of h is h∈[h_ min, h_ max] where
    
        h_ min=1/η(1+K^2_η),   h_ max=1/η .
    
    
        K_η≡W/2η .
    
    And the value of the PDF on the lower boundary is
    
        p_η(h_ min)=η^2 (1+K^2_η)^2/K_η .
    
    
    In the case of delocalized non-ergodic states  (<ref>)
    
        M≫ K_η≫ 1 .
    
    
    
    The PDF p_η(h;0)≡ p_η(h) is plotted in Fig. <ref>. 
     The PDF reaches the local  maximum on the lower  boundary h_ min corresponding to  values of
    marked state energies ≃ W located at the edges   of the IB. In the region h∼ 1/η the probability density reaches very small values, p_η(h,z)∼η^2, corresponding to the energies of marked  states  |-z|≃η. Maximum value of h=1/η  corresponds to exact resonance =z. The PDF    p_η(h;0) has an integrable singularity at this point.
    
    
    It is of interest to consider the PDF of the sum of random variables h_m over all marked states 
    
        s^h_M=1/M∑_m=1^Mh_m,   h_m=η/(z-_m)^2+η^2 .
    
    In the non ergodic  phase W≫η mean value of h_m is much smaller than its standard deviation
    
        ⟨ h_m⟩=∑_σ=± 1 arccot(2η/W-2 σ z)/W≃π/W ,
    
        ⟨ h^2_m⟩≃π/2Wη≫⟨ h_m⟩^2 .
    
    Note that the mean is dominated by small marked state energies _m∼η  while standard deviation is dominated by _m∼ W.
    
    However for sufficiently large M the mean value of the  sum ∑_m=1^Mh_m  is much greater than its standard deviation provided that ≪η
        ⟨ s^h_M⟩^2-⟨ s^h_M⟩^2/⟨ s^h_M⟩^2≃1/2π/η≪ 1
    
    Therefore in the  delocalized phase
    
        η≫=W/M ,
    
    the sum ∑_m=1^Mh_m is self-averaging.
    
    It is convenient to  introduce  rescaled variables
    
        y_m=√(h_m η) .
    
    Their  PDF has the form
    
        _η(y)=1/K_η y^2√(1-y^2) .
    
    Boundaries of the domain of _η(y) are
    
        y_ min=1/√(1+K^2_η)≤ y <y_ max=1 .
    
    
    § PDF  OF  THE  IMAGINARY PART OF SELF-ENERGY IN SELF-CONSISTENT BORN APPROXIMATION
     
    
    
    In this section  we provide details of calculations of self-consistent Born approximation presented in Sec. <ref> of the main text.
    We study the  PDF of the sum
    
        =^2∑_m=1^Mw_m η/(z-_m)^2+η^2 ,
    
     where  w_m=^2(d_0m)/^2  (see Eqs. (<ref>), (<ref>)) are   random variables sampled from the distribution g_∞(w) (<ref>) and marked state energies _m obey the dstribution p_A(/W)/W (<ref>).
     The sum in (<ref>) can be written in the form  
     
        =^2/η∑_m=1^Mx_m,   x_m=w_m y_m^2 ,
    
     where y_m are   random variables  (<ref>)  sampled from  the distribution  _η(y) (<ref>).
    For |z|≪ W random variables x_m obey  the PDF
     g_η(x) such that
    
        g_η(x)=∫_y_ min^1dy ∫_1^∞dw  _η(y) g_∞(w) δ(x-w y^2) .
    
    Using (<ref>) and (<ref>) one can show that (cf. also (<ref>))
    
        lim_η→∞g_η(x)=g_∞(x) .
    
    
    In order to calculate the  PDF of the sum  (<ref>)
     in the limit M→∞ we use GCLT following the same approach as that in Sec. <ref>.
    The PDF of the random variable  equals
    
        ()≃1/2π∫_-∞^∞dk  e^-i k +Mϕ_ηk ^2/η ,
    
    where ϕ_η(u) is the characteristic function of the PDF g_η(x) (<ref>)   
    
        ϕ_η(u)=∫_1/1+K^2_η^∞dx  g_η(x)(e^i u x-1) .
    
    
     §.§ PDF  of  individual terms in the sum
    
    
    
    After some transformations we get from Eq. (<ref>) 
    
        g_η(x)=1/x^2 K_η√(2π)
           ×∫_1/√(1+K^2_η)^min(1,√(x))dy/√((1-y^2)(log x^1/2-log y)) .
    
    The PDF is plotted in Fig. <ref>. Its maximum lies very close to the left boundary of its domain
     x∈[1/(1+K^2_η),∞ ). For x≪1 the PDF g_η(x) depends on x in terms of the rescaled parameter z=x (1+K^2_η) whose PDF is 
     
        g̅_η(z) ≃ erf√(1/2log z)/z^3/2√(2) .
    
    The plot of  g̅_η(z) is given  in Fig. <ref>), its maximum z_ max≃ 1.35. Typical values of x_m≃ z_ max/K^2_η≪ 1 correspond to w_m∼ 1 and  to a broad  PDF of marked state energies,  |z-_m|∼ W.
     
     
     
    We are interested in the limits (cf. (<ref>))
    
        x≫1,   K_η≫ 1 .
    
    We note that  log x≫ |log y| in the denominator of (<ref>) for all y except for the small interval 
    
        1/√(1+K^2_η)≤  y   ≲1/x ,
    
     whose contribution to the integral  neglected.
    Expanding the integrand in powers of (log x)^-1/2 we get
    
        g_η(x)≃π/2K_ηg_∞(x)-πlog 2/2K_ηπ^1/2x^2log^3/2x ,
    
    where function g_∞(x) is defined in (<ref>).
    
    
    
    
    We observe from  (<ref>), (<ref>) that 
    
    η⟨  h_m⟩=  π/2K_η. Using the expressions for  g_∞ (<ref>) and ⟨  h⟩ (<ref>) we obtain under the condition (<ref>) 
    
        g_η(x)≃ηh g_∞x/ηh,   x≫ 1 .
    
    Given a  large  deviation of   x_m  satisfying (<ref>),  the conditional PDF of η h_m  is  narrowly peaked around its  mean value corresponding  to  |_m-z|∼η.
    In contrast, typical values of x_m correspond to  a much broader PDF of   _m∼ W.  This gives rise to a small factor π/2K_η∼η/W in the leading order term in (<ref>).
    
    
    
    
    
    
    
     §.§ Characteristic function of the PDF of the elements in the sum
    
    
    
    
    The relation between the characteristic functions ϕ_η(u) and ϕ_∞(u) (<ref>) in the limit 
    
        |u| ≪ 1 ,
    
     should be the same as the relation (<ref>) between the corresponding 
    PDFs g_η(x) and g_∞(x) in the limit of large x (<ref>).  Here we will establish this directly.
     We break ϕ_η(u) in two parts
    
        ϕ_η(u)=ϕ_η^1(u)+ϕ_η^2(u) ,
    
    where
    
        ϕ_η^1(u)=∫_1/1+K^2_η^1dw  g_η(w)(e^i u w-1) ,
    
    
        ϕ_η^2(u)=∫_1^∞dw  g_η(w)(e^i u w-1) .
    
    Expanding  ϕ_η^1(u) in u we get
    
        ϕ_η^1(u)≃π/2K_η  i ζ_1 u ,
    
     where
     
        ζ_1=2/π^3/2∫_0^1dx/x∫_0^√(x)dy1/√((1-y^2)log(x/y^2)) .
    
    To calculate ϕ_η^2(u) in the limit of small |u| we introduce X≫1 such that
    
        |u|≪ X |u|≪ 1 ,
    
    and write
    
        ϕ_η^2(u) =ϕ_η^2,-(u) +ϕ_η^2,+(u) .
    
    Here
    
        ϕ_η^2,-(u)    =∫_1^Xdx  g_η(x)(e^i u x-1) ,
        ϕ_η^2,+(u)    =∫_X^∞dx  g_η(x)(e^i u x-1) .
    
    We use (<ref>) and expand ϕ_η^2,-(u) in u
        ϕ_η^2,-(u) ≃ i u ∫_1^Xg_η(x)dx .
    
    To calculate the term ϕ_K^2,+(u) we use the approximation (<ref>)  and write
    
        ϕ_η^2,+(u)    =π/2K_ηϕ_∞(u)
        -i u π/2K_η∫_1^Xg_∞(x)x dx .
           -πlog2/2K_η∫_X^∞dx e^i u x-1/√(π) x^2 (log x)^3/2
    
    where the characteristic function ϕ_∞ is defined in   (<ref>).
    
    Combining ϕ_η^2,±(u) together and taking the limit X→∞ we get after some transdormations  
    
        ϕ_η^2(u)≃π/2K_η(ϕ_∞(u)-i ζ_2 u)
    
        ζ_2=(32/π^3)^1/2∫_0^1dy(log(1/y)/1-y^2)^1/2 .
    
    After some transformations one can show that ζ_1=ζ_2. Therefore terms ∼ u in ϕ_η^1(u) and ϕ_η^2(u) cancel each other. Combining these two quantities together  in (<ref>) we finally get 
    
        ϕ_η(u)≃π/2K_ηϕ_∞(u)+ O( |u|/  K_η  |log u |^3/2) .
    
    As expected, this  relation corresponds to the relation (<ref>) between the PDFs g_k and g_∞.
    
    
    
     §.§ GCLT for the sum 
    
    
    
    We now  revisit the expression  (<ref>) for the PDF of the variable (<ref>)). In the limit M→∞ the
    integral over k in the r.h.s of (<ref>) is dominated by small values of the argument in ϕ_η(k ^2/η). Then using (<ref>) and (<ref>) we get after the change of a variable of integration in (<ref>)
      
        ()=1/2π∫_-∞^∞dq e^-i q/+Ω_ηϕ_∞(q/Ω_η) ,
    
    where   (<ref>) is  the characteristic  value of imaginary part of self-energy of marked states   obtained in FGR-based calculation in Sec. <ref> and quantity Ω_η equals
    
        Ω_η=π M/2 K_η=πη/ .
    
     It has a meaning of the typical number of marked states within the non-ergodic miniband of the width η
     (cf. Eq. (<ref>) and Fig. <ref>). 
     
     
     We make  a self-consistent assumption (cf. Eq. (<ref>) in the main text) and set 
    
        η= .
    
    Then, one can immediately see that
    
        Ω_η=Ω_=Ω ,
    
    where Ω is  the typical number of marked states in a mini-band  defined in (<ref>).
     
     
     Comparing the expression (<ref>) with (<ref>) and (<ref>) we represent the random variable  in the form
    
        σ_Ω x+b_Ω,  (x)=L^1,1_1(x) .
    
    Here random variable  x obeys a   Levy alpha-stable distribution    (<ref>) shown  in Fig. <ref>. The  quantities  b_Ω, σ_Ω are given below
    
        σ_Ω=√(π/4logΩ) ,
    
    
        b_Ω≃σ_Ω^-1-2/πσ_Ωlog(σ_Ω^-1) +2/π (1-γ_ Euler) σ_Ω ,
    
    Their dependence on  Ω is given in the main text, Eqs. (<ref>),(<ref>), where we should replace  M with Ω.
    
    
    
    §  NUMERICAL SIMULATIONS
    
    
    In this Section we provide details of the numerical analysis of the ensemble of Hamiltonians introduced in Sec. <ref> in addition to the results in Sec. <ref>.
    
    
    
     §.§  Numerical justification of cavity equations
    
    
    
    
    Application of cavity method to the case of the ensemble of dense matrices considered in this paper, see Sec. <ref>, exploits the similarity between the local structure of the adjacency graph of the Hamiltonian  and the Bethe lattice. The derivation of the cavity equations (<ref>),(<ref>) for the case of   outlined in Sec. <ref> neglects off diagonal terms Y in comparison to diagonal X, which is justified for graphs with extensive number of neighbors <cit.>, where,
    
        X=1/M∑ _j _1 j^2   G_jj(z), 
        
        Y=2/M∑ _j≠ k_1 j_1 k G_jk(z ),
    
    where G_ij is the single particle Green function corresponding to the Hamiltonian  at energy near the center of the band, introduced in Sec. <ref>.
    It has been shown for Levy matrices <cit.> that the ratio |Y/X| scales to zero with growing matrix size M and therefore can be neglected. This argument could be extended to PBLMs considered in this paper. We confirm the validity of this approximation numerically by analyzing the probability distribution of the ratio |Y/X| as a function of the matrix size M. In Figs. <ref>,<ref>, <ref> the distribution of |Y/X| scales towards high weight at vanishing values of |Y/X| with growing M. 
    
     
    
    
     §.§  Numerical analysis of population transfer time
    
    
      §.§.§  Population transfer time from the dynamical correlator 
    
    
    In addition to Fig. <ref> in Sec. <ref> of the main text, we perform a similar collapse of the dynamical correlator frequency dependence for different matrix sizes M for a range of different values of γ. In Figs. <ref>-<ref>  the characteristic energy scale extracted from each set of plots using this procedure Γ_ε = Γ_typ M^ε corresponds to the typical mini-band width with the respective value of the parameter γ. The fitting parameter in the scaling exponent ε is small for all γ we considered and is consistent with finite size effect.
    
    
    
    
      §.§.§   Population transfer probability as a function of time 
    
    
    In the main text we analyzed the complexity of the PT protocol using the solution of the full non-linear cavity equations for the size of the typical mini-band and estimated the number of states in the mini-band using the classical value of the level spacing W/M. In this section we analyze the scaling of the population transfer time using exact  numerical time evolution.  We contrast the population transfer time obtained from the characteristic energy scale of the frequency dependence of the dynamical correlator in Figs. <ref>-<ref> with the time dependence of the transfer probability,
    
        p(t)=| ⟨ i | ψ(t) ⟩|^2,
    
    where |i⟩ is the initial bitstring and |ψ(t)⟩ is the wave function resulting from the evolution with the impurity band Hamiltonian in transverse field , see Sec. <ref>, for a time t, which is the quantity directly observed experimentally.  Note that in Fig. <ref> the time scale at which the transfer probability becomes of order one depends strongly on the parameter γ, reflecting the fact that the characteristic time is determined by the size of the many-body mini-band Γ rather than the typical off-diagonal matrix element V_typ.  To verify this we rescaled the unit  of time with the square root of the number of states in the mini-band √(Ω), a good approximation for the scaling of the mini-band, see Sec. <ref> for qualitative discussion and Sec. <ref> for rigorous results. We observe approximate collapse of the curves for different values of γ corroborating the PT runtime scaling presented in the main text as well as the estimate of the number of states in the mini-band.
