Area law in the exact solution of many-body localized systems
Evgeny Mozgunov











 mvjenia@gmail.com
 Institute for Quantum Information and Matter,
  Department of Physics,
  California Institute of Technology









             


Many-body localization was proven under realistic assumptions by constructing a quasi-local unitary rotation that diagonalizes the Hamiltonian (Imbrie, 2016). A natural generalization is to consider all unitaries that have a similar structure. We bound entanglement for states generated by such unitaries, thus providing an independent proof of area law in eigenstates of many-body localized systems. An error of approximating the unitary by a finite-depth local circuit is obtained. We connect the defined family of unitaries to other results about many-body localization (Kim et al, 2014), in particular Lieb-Robinson bound. Finally we argue that any Hamiltonian can be diagonalized by such a unitary, given it has a slow enough logarithmic lightcone in its Lieb-Robinson bound.













                             

                              
Area law in the exact solution of many-body localized systems
    Evgeny Mozgunov
    December 30, 2023
=============================================================






§ INTRODUCTION


Physical models change their behavior in the presence of disorder  — a randomness in the terms of the corresponding Hamiltonian that's breaking the translation symmetry. First step in understanding the physics of disorder was made by Anderson<cit.>. He mathematically demonstrated a phenomenon of localization for non-interacting particles in a disordered potential. A quantum particle doesn't leave the vicinity of its starting point. This behavior of non-interacting particles can be thought of as a zeroth order approximation in the description of the interacting system. The question of convergence of any expansion around the non-interacting solution turned out to be a hard nut to crack. It was noted<cit.> that an interacting system can be thought of as a single particle traveling in the Fock space, but the corresponding potential was not localizing in the Anderson sense. First positive answer was presented by Basko et al<cit.>, which is when this problem became widely known as many-body localization (MBL). The achievement of this pioneering work was to find a realistic system (electrons in grains with a weak tunneling between them) and a range of parameters (weak interaction) where the diagrammatic expansion becomes tractable. The authors demonstrated zero conductivity at a finite temperature  — the system behaved as an Anderson insulator would.

An alternative route of investigation concentrated around disordered spin systems and conjectures by Huse et al<cit.>.

Spins-1/2 are described by Pauli matrices σ_i^α. These degrees of freedon are normally fixed, and excitations do not always admit a quasiparticle description; what localization meant was open to debate. Huse et al conjectured that dressed ("tailed" here and below) spin operators τ^α_i can be found that are conserved quantities of the dynamics. The set of those quasi-local conserved quantities τ^α_i was conjectured to form an alternative (logical) operator basis. The coefficients of the mapping between the two thus is given by the trace norm tr(σ_i^ατ^β_j) which decays exponentially according to quasi-locality of τ_i^α. This mapping allows to rewrite the spin Hamiltonian as a classical Hamiltonian (only containing τ_i^z) in the logical basis. This logical Hamiltonian contains a non-local exponentially decaying interaction inherited from quasi-locality of τ_i^α.

All of the above conjectures turn out to be true for a specific spin system solved by Imbrie<cit.>. His construction directly obtained the logical Hamiltonian as expected by Huse et al<cit.> by interleaving steps of perturbation theory with spatially sparse non-perturbative rotations. The exponential tails are adjusted as no decay happens in the regions containing those non-perturbative rotations, called resonances. But resonances are proven to be rare thus exponential decay holds on average. We will present the construction in more detail below, as it inspired our own definition of many-body localization. The path to showing quasi-locality of logical operators lies through a property called telescopic sum that we show. To our knowledge Imbrie's work<cit.> is the only first principles proof for spin systems that serves as a beacon guiding other mathematical work<cit.> in the field.

Quite a lot of progress has been made by taking a different starting point than first principles (the Hamiltonian). Most notably, a Lieb-Robinson bound has been proven<cit.> from conjectures alike ones by Huse et al<cit.>. Lieb-Robinson bound characterizes how an operator evolves under the dynamics generated by the Hamiltonian. Normally the region on which the majority of the operator weight can be found (in terms of trace coefficients with the operator basis as shown above) grows linarly with evolution time. In MBL systems the operator spreads slower, only logarithmic in time. We show that is also true for the first principles solution.

The intuition about MBL was often informed by numerical simulations, in particular Nayak et al<cit.> conjectured that a finite-depth local circuit can map an eigenstate approxiamtely to a product state. We will confirm that conjecture by providing a bound on the error of such mapping, and emphasize that a single unitary works equally well for all eigenstates. A circuit on top of a product state can be thought of as a matrix-product state. Some of the authors focused on a so called strong MBL, where they were able to prove<cit.> a matrix-product state representation of every eigenstate. They also conjectured that a stronger Lieb-Robinson bound where the lightcone does not spread at all holds for the states at low energies. We find no evidence for that claim. 









§ PROPERTIES OF THE FAMILY OF IMBRIE'S SOLUTIONS

Every Hamiltonian H is characterized by its eigenvalues H_diag and eigenvectors U:

    UHU^†=H_diag

 Imbrie<cit.> presented a construction of U = ∏_k U_k that works for a specific Hamiltonian on a 1d chain:
 
    H=γ∑_iσ^x_i +∑_ih_i σ^z_i + J_i σ^z_i σ^z_i+1

 with -1 ≤ h_i,J_i≤ 1 uniform i.i.d. random numbers and γ≪ 1. The structure of  U = ∏_k U_k is illustrated in Figure <ref>. We shall present the construction abstractly, without specifying all of the internal structure of U_k. We believe that for any MBL system a generalization of Imbrie's construction can be defined, so our presentation can be read both as a restatement of Imbrie's proof, and as a postulate about any system sufficiently deep in the MBL phase.


First consider no resonances (no red blocks in Figure <ref>). There are infinitely many terms in the product  U = ∏_k U_k counted by k=1,2... 
(Support)
Each block
U_k =e^A^(k),
 where iA^(k) is a Hermitean generator
 that is a sum of floor(15/7L_k)-local terms on each site of the chain, where L_k =(15/8)^k:

    A^(k) = ∑_i A^(k)_i


Here locality is in terms of the interval on the chain that operator covers  — the distance in space between the leftmost and the rightmost spin that the operator depends on. We provide a guide to the text of Imbrie's paper in the Appendix B, explaining how we chose the number 15/7 in the support definition above. We also conclude that Imbrie's proof allows one to restrict the norm of A^(k)_i as follows:
 (Norm)

    A^(k)_i≤ c'χ^ceilL_k-1

where χ = c γ^0.95 and c, c' are constant numbers.

Note that the magnetic field γ enters the definition here. The construction rests on the fact that the 15/7L_k appearing in the support has the same scaling with k as L_k-1 = 8/15L_k. The coefficient c is obtained as a polynomial of exp(suppA^(k)_i/L_k). If scalings were to be different even by a logarithmic factor (say, supports grow as kL_k), then instead of (cγ^0.95)^L_k-1 the bound would contain L_k-1! γ^L_k-1, which diverges regardless how small γ is. In Imbrie's paper the fact that support has ∼ L_k scaling is used several times, however the proof of it raises some questions, as we note in Appendix B. These two postulates are enough to prove the convergence of the no-resonance infinite product  U = ∏_k U_k in the following sense:

Operators acquire exponentially decaying tails when conjugated by no-resonance U.

    UXU^†-U_cXU_c^†≤ 24· 2^0.5S(X)α^c+1X,
    U^†XU-U_c^†XU_c≤ 24· 2^0.5S(X)α^c+1X

where U_c is constructed using the generators A^(k) truncated to a collar of c sites of the support of X (all the terms reaching outside the collar are dropped). α = 24c'χ^ 56/225 is a power of transverse field γ.

For the proof, see Appendix A.

    
    < g r a p h i c s >

    
    < g r a p h i c s >


Consider the implications of this result for the eigenvalues H_diag. 

    H_diag =UHU^†=γ∑_i(σ^x_i)_tailed +
    +∑_ih_i (σ^z_i)_tailed + J_i (σ^z_i σ^z_i+1)_tailed

On the other hand, since H_diag is a diagonal matrix:

    H_diag =f(σ_z) = ∑_s⊂{1… L} J_s ∏_i∈ sσ_i^z.

The bound on J_s can be directly inferred from the general bound on the telescopic sum (Theorem <ref> in the Appendix):


    J_s ≤ 12 α^d(s)

where α = 24c'χ^ 56/225,  d(s) is consecutive support and no resonances are within it.


Using a time-averaged version of this (see Section <ref>) as a starting point in an ensemble-averaged sense, Kim et al<cit.> derive Lieb-Robinson bound:


    [A,B(t)]≤ ate^-bd_A,B=
     = atα^d_A,B

where B(t) = e^iHt B e^-iHt, α = 24c'χ^ 56/225, a=cAB and c is a constant number.

We provide an independent derivation in Appendix A. Another known MBL effect of logarithmic growth of entanglement across the cut can also be proven<cit.>. We prove that U is close to a finite-depth local circuit in Appendix <ref>.

The benefit of taking Postulate <ref> as a starting point is that one gets area law essentially for free. Consider an eigenstate |E_s⟩= U^†|s⟩  where |s⟩ is a product state of spins in z-basis and a cut at some location:


    S_cut = S(ρ_left)< 650 c'χ

where ρ_left = tr_right |E_s⟩⟨ E_s|.

No other proofs of area law in an eigenstate appeared in the literature to our knowledge. It was considered to be an obvious corollary of exponential decay of correlations, however Imbrie did not prove exactly the same decay of correlations as required for 1d area law.<cit.> We show our proof in Appendix A, and obtain the Postulate <ref> from the construction by Imbrie<cit.> in Appendix B.

We now interleave the circuit with resonances  — locations where perturbation theory in γ failed and a local non-perturbative rotation had to be applied.
 (R-size)
At each step, the operators acting between e^A^(k) are 

    O_k =∏_r O_k,r,

where O_k,r are big unitary rotations supported on <4.2L_k =4.2(15/8)^k  consecutive spins.

Let R_k be the set of spins involved in O_k. Imbrie<cit.> explicitly presents a following bound on how rare they are:
 (R-density)
The density of step-k resonances

    P(i∈ R_k) = lim_L→∞|R_k|/L≤ cϵ^c'log^2L_k

where ϵ = γ^0.05 and c, c' are constant numbers.

We postulate 1-4 as a definition of a full U. Now entanglement bound strongly depends on the position of the cut!


    S_cut≤   650 c'χ + 4ϵln2

and max(S_cut)∼log L 

Here S_cut is the average over different positions in the system (equivalently, different realizations of disorder), and the true finite size scaling is more complicated than logL but is slower than any power. We note an incredibly small power in the resonance density bound ϵ = γ^0.05. That means that Imbrie's construction as is cannot be applied for magnetic fields γ≥ 10^-6. MBL survives for much bigger magnetic fields, and this discrepancy is due to multiple overestimations needed for Imbrie's proof to got through. If one uses a more strict choice of which regions to call resonant, their density remains controlled roughly up to the transition to extended phase. We demonstrate it numerically in a companion paper<cit.>.



§ TIME-AVERAGING


Consider the following transformation of a local operator X:

    X' =lim_T→∞1/T ∫_0^T e^iHt X e^-iHt dt;

this kills off-diagonal elements in the eigenbasis of H, so commutation with H is enforced.
Time-averaging is used in Kim et al<cit.> as part of a definition of MBL. They show the following for time-averaging with an MBL Hamiltonian:

Infinite time-averaging introduces exponential tails.

    X'-X'_c≤ a(4α)^c

where X'_c is the truncation of X' to a region with collar c around the support of X, α = 24c'χ^ 56/225 and a,Z are constant numbers.

 
Following<cit.>, logical spin operators Uσ^α_i U^† form an operator basis and obey the telescopic sum bound proven in Appendix A. Coefficients of an operator X in this basis are given by traces of products trXUσ^α_i U_†, which are bounded as c'α^i-i_X where i_X is the nearest site in X. Only the coefficients of combinations of Uσ^z_i U^† survive the infinite time averaging. Then each term has a telescopic bound c'α^c where c is the size of the collar. Two bounds combine to give a collar bound on X'. The α is only increased by a square of local dimension as there are many possible operator strings.
 


If we time-average a local patch of the Hamiltonian H_a in H= H_a +H_b  for a finite-time, the commutation with H is not fully enforced. Let's bound the commutator of the two parts:

    [H_a',H_b'] =[H_a',H] = 1/T∫_0^T e^iHt [H_a,H] e^-iHt dt  =
    
        =1/iT∫_0^T d/dte^iHt H_a e^-iHt dt  =e^iHT H_a e^-iHT - H_a/iT
     [H_a',H_b']≤2S/T,

where S is the support of H_a. The spread of the local operator of any Hamiltonian can be bounded by the Lieb-Robinson bound. Let the light cone (where the bound on [A(t),B]∼ 1) be given by x(t), therefore time-averaging for time T should give an operator whose support outside the radius x(T)+Δ x can be bounded as e^-cΔ x, where c is a Hamiltonian-independent number for 2-local Hamiltonians (and can be strengthened even more for MBL Lieb-Robinson bound: e^-cΔ x = α^Δ x now H-dependant unlike the generic Lieb-Robinson bound). Moreover, the time-averaging can be done with the Hamiltonian truncated to an x(T)+Δ x collar of the local operator, with only exponentially small error:

    H̃_a' =1/T ∫_0^T e^iH_ct H_a e^-iH_ct dt =   H_a' + δ_a,   δ_a≤ e^-cΔ x.

 Note that δ is a nonlocal interaction that's exponentially small in the distance as e^-cd. Now split the system in many regions a_1,a_2 … So if we take H_a to have support S>x(T)+Δ x, the effects of the time-averaging could be truncated in such a way that H̃_a_i' only overlap with their neighbors H̃_a_i±1'. So we have reduced our Hamiltonian to a nearest neighbor chain of approximately commuting terms (as in Bound <ref>), with an exponentially small error:
 
    H = ∑_i H̃_a_i' + δ_i.

Now let's take even-odd pairs H̃_a_2k'+ H̃_a_2k+1'. We would like to use a result by <cit.> that reduces approximately commuting operators to exactly commuting ones, by small adjustments to each:

    H̃_a_2k' =A_2k+ϵ_2k,   H̃_a_2k+1' =A_2k+1+ϵ_2k+1, 
    
         [A_2k,A_2k+1] =0,

and the ϵ<(2S/T)^1/6 and have the same support as the overlap of a_2k' and a_2k+1' (the supports of corresponding H'). Now we introduce the notation

    h_k = A_2k + A_2k-1,     [h_k,h_k+1]=0, 
       H= ∑_k h_k +∑_i (ϵ_i + δ_i).

We have reduced any local Hamiltonian on a chain to a nonlocal one with exponentially decaying long-range interaction, written in the form of commuting terms plus a small perturbation. We note that to repeat this in 2d one would need a version of Hastings' result for three matrices, and that one has a counterexample presented in his paper. 

For chains, one may think of applying Imbrie's construction to this form of the Hamiltonian - after step 1, Imbrie's Hamiltonian has exactly this form with long range but rapidly decaying interactions! Unfortunately, the bounds for a general 1d Hamiltonian do not work out in our favor: the number of levels in the block S>x(T)+Δ x grows as 2^S, so the minimum level distances decrease as 2^-S. However, the norm bound on ϵ is just a power 1/T. For the standard Lieb-Robinson bound x(T) =vT and S>x(T), so the level distances are always smaller than the perturbation scale, which is the opposite to what is used in Imbrie's construction. However, in MBL with x(T) ∼lnT, the scales may just work out: T ∼(1/α)^x(T), so the norm of the perturbation goes as α^c̃S where S is our choice of the support and c̃ is some constant. It decreases faster than (1/2)^S for a sufficiently small α. We have shown the following:












If H obeys MBL Lieb-Robinson bound

    [A,B(t)]≤ atα^d_A,B,

then the perturbation on top of a S-local commuting component of H
 V=∑_i V_i is sufficiently small for Imbrie's proof to go through:

    V_i≤ c”α^c'S

where c”,c', S are constant numbers.

In this way, any Hamiltonian that possesses the MBL Lieb-Robinson bound can be reduced to a form "diagonal + small perturbation", even if the z-basis is not apparent from the start.
*




§ PROOFS OF PROPERTIES OF IMBRIE CIRCUIT



 §.§ Telescopic sum for Imbrie circuit



  §.§.§ Small unitary expansion
 We consider a generator A^(k)  of a unitary evolution e^A^(k) (corresponding to step k) such that A^(k)= ∑_i A_i^(k) where A_i^(k)≤ c' χ^ceilL_k-1) and each A_i is supported on floor(8/7L_k+1) consecutive sites. We would like to study the structure of a local operator X, evolved as e^A^(k) X e^-A^(k). For the purposes of the proofs below, we use the term support S(X) of any operator X as the shortest interval of the chain containing all the sites where action of X is distinct from identity. We prove




    e^A^(k) X e^-A^(k)  = ∑_j=0^∞ X_j

where  S(X_j) ≤ S(X) +  jΔ S,  Δ S =2(floor8/7 L_k+1-1)  and 

    X_j ≤ 2^j_0(2c'Δ S )^j χ^jceilL_k-1X

where  j_0 = ceilS(X) +S(A) -2/Δ S and S(A) = floor8/7 L_k+1




So the operator X acquires exponential tails after conjugating with e^A^(k).



We drop the superscript (k) in A^(k) for convenience.

    e^A X e^-A = ∑_n A^n/n! X ∑_m (-A)^m/m!

collecting the terms that have the same power of A

    =∑_j ∑_m=0^k A^j-m/(j-m)! X (-A)^m/m!

Compare this to the expansion of j'th order commutator

    [A[A…[A,X]]]^(j) = ∑_m=0^j C_m^j A^j-m X (-A)^m

where C_m^j = j!/(j-m)! m!

    ∑_j ∑_m=0^j A^j-m/(j-m)! X (-A)^m/m! = ∑_j 1/j![A[A…[A,X]]]^(j)

Now denote X_j =  1/j![A[A…[A,X]]]^(j). We found that

    e^A X e^-A  = ∑_j X_j


Recall that support of X is S(X) consecutive sites and the support of A is S(A) = floor8/7 L_k+1. A commutator [A,O] for a local operator O of support S(O) consecutive sites is nonzero for only the terms in A = ∑_i A_i that overlap with O. The index i goes over every site, so there are S(O) + S(A) -2 such terms. The support of [A,O] becomes S(O) + 2S(A)-2.  The support of X_j which is the j'th order of the commutator just grows linearly:

    S(X_j) =S(X) + j(2S(A)-2) = S(X) + jΔ S

where we denoted Δ S = 2S(A) -2. The number of terms of A that have nonzero commutator at j'th step is given by S(O) +S(A)-2 for O=X_j-1:

    n_j = S(X_j-1) + S(A) -2 = S(X) +jΔ S +S(A) -2

For the norm bound on X_j, a prefactor before the j'th power of A is 1/j! and a product of all n_j's:

    X_j≤1/j!(∏_i=1^j n_i )A^j X

Let us bound the prefactor:

    1/j!∏_i=1^j n_i = 1/j!∏_i=1^j (S(X) +S(A) -2+iΔ S ) =Δ S^j/j!∏_i=1^j (S(X) +S(A) -2/Δ S+i ) ≤
    ≤Δ S^j/j!∏_i=1^j (ceilS(X) +S(A) -2/Δ S+i ) =  Δ S^j/j!∏_i=1^j (j_0+i ) = Δ S^j (j+j_0)!/j!j_0! = Δ S^j C_j_0^j+j_0

where we denoted j_0 = ceilS(X) +S(A) -2/Δ S. Finally we use C_j_0^j+j_0≤ 2^j+j_0 to arrive at

    1/j!∏_i=1^j n_i ≤ (2Δ S)^j 2^j_0

Combining that with the bound on A, we get the desired form:

    X_j ≤ 2^j_0(2c'Δ S )^j χ^jceilL_k-1X







  
Aside: comparison with Lieb-Robinson Bound
We have found that if one conjugates a local operator X with S(A)-locally generated 1d evolution e^A for small time t≤ c'χ^ceilL_k-1, operator X acquires exponential tails.


It is a much stronger result than if we apply Lieb-Robinson Bound for small χ for propagation of tails beyond the collar, but becomes trivial for χ=1. Here's this LRB estimate for comparison - it still gives nontrivial answers even for big χ:
 
    ∑_j=c^∞X_j ≤exp{-(ac -  S(A)c'χ^ceilL_k-1}X

Here a is just some number <20 much smaller than ceilL_k-1/(-lnχ).

We shall use the stronger first bound as in our case χ≪ 1.



  §.§.§ Infinite product of unitaries
 Now we'd like to use our bound to estimate precision of collar approximations to:

    ∏_k=1 e^A^(k) X ∏_k'=1 e^-A^(k')

 as in the total Imbrie circuit without resonances. The order in the product will turn out to be unimportant for our bound. From the point of view of the physical meaning, if smaller k act on the operator X first, then X is a physical operator and evolved X represents the action on "logical" spins, or integrals of motion. In particular, for the diagonalization of the Hamiltonian, smaller k act on the physical Hamiltonian first. If smaller k act last, then X is a logical operator, and its evolved form will be the representation of such operator on physical spins. 
 
 
 We now present and prove the bound:
 



    ∏_k=1^∞ e^A^(k) X ∏_k'=1^∞ e^-A^(k') = ∑_j=0^∞ X_j,    S(X_j) =S(X) +2j,   Y_n≤ 6· 2^0.5S(X)α^j X

where α = 24c'χ^ 56/225 and the order in the product is unimportant.

 If we add a collar of c sites to the support of X, and take A_c^(k) to be the sum of terms fully within the collared region, then 

    ∏_k=1^∞ e^A^(k) X ∏_k'=1^∞ e^-A^(k') - ∏_k=1^∞ e^A^(k)_c X ∏_k'=1^∞  e^-A^(k')_c≤ 24· 2^0.5S(X)α^c+1X

 α≤ 0.5 is required for the second bound to hold.


 



 
 We start from looking at individual rotation, that was discussed above:
 
    U_k X U_k^†  = ∑_j_k X_j_k,     S(X_j_k)  =  S(X) +2j_kfloor(8/7 L_k+1-1),     X_j_k≤ 2^j_0(2c'Δ S )^j χ^jceilL_k-1X

 If there are two rotations, we get:
 
    U_k' U_k X U_k^† U_k'^† = ∑_j_k,j_k' X_j_k,j_k', 
      S( X_j_k,j_k')  =  2j_kfloor(8/7 L_k+1-1) + 2j_k'floor(8/7 L_k'+1-1),  
     X_j_k,j_k'≤ 2^j_0,k + j_0,k'(j_k)(2c'Δ S_k )^j_k(2c'Δ S_k' )^j_k'χ^j_kceilL_k-1 +j_k'ceilL_k'-1X

we see that the order of unitaries indeed does not matter except for the first term, that depends on the support of previous terms. We will see in paragraph <ref> how that term will remain bounded for every X_j_k,j_k' s.t. S(X_j_k,j_k') =S(X) +2j regardless of the order. For an infinite product, the components of evolved X are indexed by a string of integers {j_k} for all k=1,2,3…: 

    ∏_k=1 e^A_k X ∏_k'=1 e^-A_k' = ∑_{j_k} X_{j_k}, 
      S(X_{j_k})  = S(X)+  ∑_k2j_kfloor(8/7 L_k+1-1) ,  
     X_{j_k}≤ C({j_k},order) (∏_k (2c'Δ S_k )^j_k) χ^∑_kj_kceilL_k-1X
    
    C({j_k},order) =2^j_0,k(1) + j_0,k(2)(j_k(1)) +…

Terms X_{j_k} whose resulting supports S( X_{j_k})=S(X) +2j (j∈Z)  are combined into X_j:

    X_j = ∑_{j_k}|S( X_{j_k})=S(X) +2jX_{j_k}

We will provide a uniform bound in terms of j on the norm of each X_{j_k} in X_j, and also argue that the total number of terms is only exponential in j.

One of the terms contributing toX_{m_k}|j is  

    A =∏_k(2c'Δ S_k )^j_k

We know that Δ S_k = 2floor(8/7 L_k+1-1). We will leave one of the powers of two to the C({j_k},order), and bound a 

    A =∏_k(2c'floor(8/7 L_k+1-1))^j_k≤∏_k(16/7c' L_k+1)^j_k

The collar is increased in steps of S(A_k)-1 =floor8/7 L_k+1 -1 on each side, so 

















    j= ∑_k (S(A_k)-1)j_k

We can collect two sums over j_k in our bound:

    A ≤(16/7c')^∑_k j_k(15/8)^∑_k (k+1)j_k

Note that

    ∑_k j_k ≤∑_k (k+1)j_k ≤∑_k (S(A_k)-1)j_k = j

so assuming c'>1 (or else it can be replaced by 1 in the original bound):

    A ≤(30/7c' )^j

That's a good enough bound for our purposes. Now let's turn our attention to powers of χ:


    χ^∑_kj_kceilL_k-1

We'd like to bound it by a j-dependent expression as well. χ<1, so we need to bound the sums in the opposite direction:

    ∑_kceilL_k-1j_k ≥∑_kL_k-1j_k  = (8/15)^27/8∑_k8/7L_k+1j_k ≥
    ≥(8/15)^27/8∑_kfloor(8/7L_k+1-1)j_k = (8/15)^27/8 j

So the resulting power of χ is at least

    χ^∑_kj_kceilL_k-1≤χ^ (8/15)^27/8 j

Finally, we deal with the order-dependent factor:



  
2^j+j_0 factor
The factor 2^j_0 might be problematic for the proofs involving the full circuit, as X is oftentimes a big object, and the earlier steps expand it even bigger. 


For the proof above, we needed the telescopic sum bound:

    → X_j_k≤ 2^j_k+j_0,k(c'Δ S )^j_kχ^j_kceilL_k-1X

For convenience, we pulled out 2^j. All factors except the first are independent on the multiplication order and their contributions were bounded above. We will show that the product of 2^j_k+j_0,k over all steps is bounded by M^j for some constant M. The bound is not multiplication order-dependent.

The contribution of terms 2^j_k+j_0,k to a single representative of the X_j is:

    ∏_k=1^k_max 2^j_k+j_0,k = 2^∑_k=1^k_maxj_k+j_0,k

 Now we need to bound the sum:

    ∑_k=1^k_maxj_k+j_0,k =  ∑_k=1^k_maxj_k+ceilS(X_k-1) + S(A_k)-2/2S(A_k) -2

where we have introduced X_k-1 as the term that we took from the previous step. The sequence of supports S(X)=S(X_0),S(X_1)… S(X_k_max)=2j + S(X_0) should be monotonic. S(A_k) =floor8/7 L_k+1. Finally, j_k is the number of terms in the commutator we selected at step k, so S(X_k) = S(X_k-1) + j_k(2S(A_k)-2), so j_k = (S(X_k) - S(X_k-1)) /(2S(A_k)-2)

    ∑_k=1^k_maxj_k+j_0,k≤∑_k=1^k_maxS(X_k) + S(A_k)-2/2S(A_k) -2 +1 ≤3/2 k_max +(S(X_k_max)-1) ∑_k=1^k_max1/2S(A_k) -2

The series at the end is bounded:

    ∑_k=1^k_max1/2S(A_k) -2  = 1/2∑_k=1^k_max1/floor8/7L_k+1 -1≤1/2∑_k=1^∞1/8/7L_k+1 -2≤1/2∑_k=1^∞1/8/7L_k = 7/301/1-8/15= 1/2

So we arrive at a bound:

    ∑_k=1^k_maxj_k+j_0,k≤  j + 3 k_max  + S(X_0) -1/2

we note that k_max∼lnj, so we indeed only get an extra power of j, which will be suppressed by a small enough χ. Let's recall the specific dependence of k_max on j. X_0=X, then X_1 is everything within a one site collar, X_j - within a j-site long collar. The collar is increased in steps of S(A_k)-1 =floor8/7 L_k+1 -1 on each side, so floor8/7 L_k_max+1 -1≤ j:

    k_max = max_k  s.t. floor8/7 L_k+1 -1≤ j
    8/7 L_k_max+1 -2≤ j,    
     (15/8)^k_max=L_k_max≤7/15(j+2)  
    
     k_maxln(15/8)≤ln (j+2) + ln(7/15)    k_max≤ln ((7/15)(j+2))/ln(15/8)≤ (7/15)(j+2)/ln(15/8)

since 7/15(j+2)>1. So our previous bound becomes:

    ∑_k=1^k_maxj_k+j_0,k≤(2 + 7/5ln(15/8)) j  + 14/5ln(15/8)-1 +S(X_0)/2≤ 2.2 j + 1.8+0.5 S(X_0)

And the prefactor is bounded as:

    ∏_k=1^k_max 2^j_k+j_0,k≤ 2^2.2 j + 1.8+0.5 S(X_0) = (2^2.2)^j 2^1.8 +0.5S(X_0)


Now putting it all together, we've learned that

    X_{j_k}≤ (2^2.2)^j 2^1.8 +0.5S(X_0)(30/7c' )^jχ^ (8/15)^27/8 jX

for X_{j_k}∈ X_j. To bound X_j = ∑_{j_k}|jX_j we need to upper bound the number of terms in the sum.



  
Number of terms

 



 


We need to estimate the number of terms with collar j that are generated in our telescopic sums ∑_{j_k} X_{j_k}. First note that terms with j_k 0 for at least one k>k_max do not contribute. We're bounding the number of ways 

    ∑_k j_k (floor8/7 L_k+1-1) =j

can be achieved. Note that for each j_k

    j_k (floor8/7 L_k+1-1) ≤ j

So there are at most j/(floor8/7 L_k+1-1) choices for each j_k up to k_max. This gives an upper bound on the number of terms:

    ∑_{j_k}|j1 ≤∏_k=1^k_maxj/floor8/7 L_k+1-1 = j^k_max/∏_k=1^k_maxfloor8/7 L_k+1-1

For the upper bound on ∑_{j_k}|j1, we need to lower bound the denominator:

    ∏_k=1^k_maxfloor8/7 L_k+1-1 ≥∏_k=1^k_max L_k = (15/8)^∑_k=1^k_maxk = (15/8)^k_max(k_max+1)/2

We obtain the bound on the number of terms as:

    ∑_{j_k}|j1 ≤ j^k_max(8/15)^k_max(k_max+1)/2

Note that we will need the opposite inequality on k_max:

    k_max = max_k  s.t. floor8/7 L_k+1 -1≤ j
    floor8/7 L_k_max+2-1 >j
    
     L_k_max+2 > 7/8 j
    (15/8)^(k_max+1) > 7/15j

So

    ∑_{j_k}|j1 ≤ j^k_max(7/15j)^-k_max/2 = (15/7j)^k_max/2

Using the inequality k_max≤ln ((7/15)(j+2))/ln(15/8) derived before, we obtain:

    ∑_{j_k}|j1 ≤exp(ln(15/7j)k_max/2) ≤exp(ln(15/7j)ln ((7/15)(j+2))/2ln(15/8))

Expression on the right-hand side goes as exp(ln^2j), so it is clearly bounded by some power of e^j. To find the specific factor, we bound

    ln∑_{j_k}|j1 ≤ln(15/7j)ln ((7/15)(j+2))/2ln(15/8) = lnj ln(j+2) + ln15/7ln(1 +2/j) -  ln^215/7/2ln(15/8)≤
    ≤ln^2(j+2)+ ln15/7ln3 -  ln^215/7/2ln(15/8)

now we use the property of a function ln^2 x≤ 0.6x for x≥ 3. As an exercise, let's prove it. First note that second derivative of ln^2x changes sign at lnx =1, x =e  and is negative for x>e, in particular for x≥ 3. Next we calculate the values of ln^2x and 0.6x as well as their derivatives at x=3, which allows us to prove that  ln^2 x≤ 0.6x at least until x=7. Then we evalute the slopes at x=7 and see that the slope of ln^2 x is 0.55<0.6. Since the second derivative is always negative, the ln^2 x≤ 0.6x will be satisfied for all x>7. A numerical check shows that the inequality can be tightened to ln^2 x≤ 0.55x, but we will use the rigorous 0.6:

    ln∑_{j_k}|j1 ≤0.6(j+2)+ ln15/7ln3 -  ln^215/7/2ln(15/8) = 0.3j+ 0.73/ln(15/8)

Thus

    ∑_{j_k}|j1 ≤ (15/8)^0.3j+ 0.73

Now we have all the ingredients for the bound on the X_j in the telescopic sum:

    X_j≤∑_{j_k}|jX_{j_k}≤max_{j_k}|jX_{j_k}∑_{j_k}|j1 ≤
    ≤  (2^2.2)^j 2^1.8 +0.5S(X_0)(30/7c' )^jχ^ (8/15)^27/8 j(15/8)^0.3j+ 0.73≤ 6· 2^0.5S(X)(24c'χ^ 56/225)^jX




  
Final step Note that with the bound telescopic sum we can derive a bound for collared X̃ =∏_k e^A_c^(k) X ∏_k'e^-A_c^(k'), where c is the extra support on each side of X that is allowed to appear in X̃. To that accord, A_c^(k) for every k is chosen to include only terms in A that lie fully within a collar of c around X. Clearly A_c^(k) obeys the same bounds as A^(k), so in particular the decomposition X̃ = ∑_j X̃_j obeys the bound (<ref>). We see that ∏_k e^A^(k) X ∏_k'e^-A^(k') and X̃ only differ in terms starting from X_c+1 (because they involve A_i^(k)'s in regions where X̃ is different). So the difference is bounded by the sum of those:

    →∏_k e^A^(k) X ∏_k'e^-A^(k') -∏_k e^A_c^(k) X ∏_k'e^-A_c^(k')≤∑_j=c+1^∞X_j + X̃_j≤ 2∑_j=c+1^∞X_j

The sum converges:

    ∑_j=c+1^∞X_j≤∑_j=c+1^∞ 6· 2^0.5S(X)(24c'χ^ 56/225)^jX =  6· 2^0.5S(X)(24c'χ^ 56/225)^c+1/1 -24c'χ^ 56/225X

for 24c'χ^ 56/225 <1 which requires very small χ≲10^-7/c'^4 due to the untight bounds used in the derivation. Note that the original Imbrie's construction also required at least χ^1/20 <0.5 which results in χ≲ 10^-6  — it's not new for such proofs to accumulate large numbers. We can bound the denominator by 2 for  24c'χ^ 56/225 <0.5:
The sum converges:

    ∑_j=c+1^∞X_j≤ 12· 2^0.5S(X)(24c'χ^ 56/225)^c+1X

which results in the bound for the collared approximation:

    ∏_k e^A^(k) X ∏_k'e^-A^(k') -∏_k e^A_c^(k) X ∏_k'e^-A_c^(k')≤  24· 2^0.5S(X)(24c'χ^ 56/225)^c+1X










 §.§ Circuit approximation




We would like to have a FDL circuit approximation of a unitary e^A. 
Consider a simple case of S(A)=2 first: A = ∑_i A_i,i+1. The Trotter approximation in N steps is defined as:

    U_T =(e^1/NA_evene^1/NA_odd)^N

where A_even = ∑_i   even A_i,i+1. The error is given by the next term in BCH formula while it remains small:

    e^A = U_T + O( [A_even,A_odd]/N)=U_T + O(L A_i,i+1^2/N)

 Note that the error over whole macroscopic system diverges with size. However, over a small patch the error is small. Let's see how the divergence arises. Here we use big O(x) notation in a sense that e^A -U_T≤ Mx for some range of x around 0. To get the actual value of constant M in the bound, one way is to expand the following combination of unitaries e^ϵ, e^δ:
 
    e^ϵ e^δ -e^ϵ + δ = [ϵ,δ]/2 + ∑_j=3^∞ϵ^j+δ^j - (ϵ+δ)^j/j! + ϵ∑_j=2^∞δ^j /j!+ δ∑_j=2^∞ϵ^j /j! + ϵ/2∑_j=1^∞δ^j /j!+ δ/2∑_j=1^∞ϵ^j /j!

 Taking the norm of the above expression, one arrives to multiple series like ∑_j=2^∞ϵ^j /j! = e^ϵ-1-ϵ. To upper bound this function one uses a version of Taylor expansion with a mean-value form of the remained:
 
    e^x = 1+ x + e^x̃x^2/2  where  x̃≤ x

With that we will be able to derive 
  
    e^ϵ e^δ -e^ϵ + δ = γ,  γ≤ M [ϵ,δ] +O(ϵ^3,δ^3, ϵ^2δ,ϵδ^2)

And the number in front of every term can be obtained explicitly. In fact, the higher order terms are also combinations of commutators (e^ϵe^δ is an element of a Lie group that can be defined free from matrix multiplication), so higher powers of L only appear as L^pA_i^2p or lower, but we will not need it for our purposes here. We now plug in e^1/NA_evene^1/NA_odd - e^1/NA = γ where γ is bounded as above into the Trotter formula:

    e^A -U_T = (1+γ)^N -1 ≤ e^γ N-1

where we have used that all the coefficients in the Taylor expansion of e^γ N are strictly bigger than ones in (1+γ)^N. Now we can again use the mean-value form of the remainder of the Taylor expansion to bound the above with explicit constant as long as e^γ N is not too big:

    e^A -U_T≤ O([A_even,A_odd]/N)  = O(L A_i,i+1^2/N)

We will stick with the big O notation for convenience. Let's prove that for an operator X supported over S(X) sites the conjugation with e^A is simulated by our Trotter approximation with system-size independent error:


for A= ∑_i A_i,i+1

    e^A X e^-A = U_T X U_T^† + O ( ||X|| (S(X)+4N)/N||A_i,i+1||^2 )


We see that the system size dependence drops out, also U_T X U_T^† is an operator supported on S(X)+4N sites, and expressed by a tensor network of size (S(X)+4N) × 4N. The terms outside the causal cone of X just cancel.

First we use the bounds like those in Lemma <ref> for e^A X e^-A:

    e^A X e^-A =∑_j X_j,    S(X_j) =S(X) +2j,   X_j≤ 2^0.5S(X) (4A_i,i+1)^j X

Then we use the approximation by collared evolution like the one in Theorem <ref> (the final part of the proof) 

    e^A X e^-A - e^A_c X e^-A_c≤ 2 ∑_j=c+1^∞X_j≤  4·2^0.5S(X) (4A_i,i+1)^c+1X

Here c is the number of sites we are to include in the "collar" of X. A_c are the terms in A supported inside the collared region.  



The problem reduced to finding an FDL approximation to e^A_c supported on S(X) + 2c - a finite number of spins. For that purpose, Trotter approximation is used. 





The total precision defined as: 

    e^A X e^-A - e^A_c_T X e^-A_c_T ≤   4·2^0.5S(X) (4A_i,i+1)^c+1X  +O((2c+S(X)) A_i,i+1^2X/N)





Note that if we were Trotter approximating e^A instead of the collared one, we still get the same result e^A_T X e^-A_T = e^A_c_T X e^-A_c_T by cancellations of the spare terms as long as c≥2 N. So there exists an FDL circuit e^A_T over the whole system that approximates e^A locally in every patch.

Plugging in the collar as small as c=2 , we find the first term ∼(4A_i,i+1)^3 to be negligible (we assume 2^0.5S(X) is suppressed by powers of A_i,i+1 for big O bound to make sense). The only allowed N for c=2 for the total circuit to reduce to collared circuit is N=1 The error is now:

    e^A X e^-A - e^A_c_T X e^-A_c_T ≤   O((4+S(X)) A_i,i+1^2X)

More generally, we can take bigger c=2N, and the first term will always be suppressed as ∼(4A_i,i+1)^2N+1, which allows to neglect it compared to the second term:

    e^A X e^-A - e^A_c_T X e^-A_c_T ≤   O((4N+S(X)) A_i,i+1^2X/N)


Note that the precision never really gets better than O(A_i,i+1^2X). That's concerning because our original collar c=2 approximation already had ∼(4A_i,i+1)^3 precision. So if we want to approximate that well, we need to resort to a higher order Trotter scheme (instead of e^1/NA_evene^1/NA_odd it uses a more complicated step). Alternatively, since the dimension of the space we act upon (d= 2^S(X) +4 for c=2) is relatively small, we can instead use the Solovay-Kitaev theorem<cit.>. For this particular case it can be phrased as follows: to approximate a unitary on a d-dimensional Hilbert space to a precision ϵ, one needs n=O(d^aln^b1/ϵ) 2-qubit gates acting on nearest neighbors in 1d. Here a,b are some constant numbers. So to achieve precision ϵ =A_i,i+1^c+1 , one would need n=O(c2^ac) gates, or O(2^ac) depth. The number of gates in the c'th order of Trotter step also scales exponentially with c, giving the same A_i,i+1^c+1 precision. So whichever method wins depends on the constant a and the scaling of the size of Trotter step.




Now let's apply this logic to the Imbrie circuit ∏_k=1^∞ e^A^(k). We will see that as early as for approximation up to k= 2 the first order Trotter scheme becomes useless. Indeed, the S(A_k)=floor(8/7L_k+1), so S(A_1) = 4, S(A_2) =7, S(A_3) =14. The telescopic sum structure X= X_j allows us to drop higher order steps, since step j only uses the terms with S(A_k)-1≤ j, and make an error given by Theorem <ref>:

    ∏_k e^A^(k) X ∏_k'e^-A^(k') -∏_k^k_max e^A_c^(k) X ∏_k'^k_maxe^-A_c^(k')≤  24· 2^0.5S(X)(24c'χ^ 56/225)^c+1X

For c<6,   S(A_2)-1=6>c, so we only include the 1st order term in the product and obtain the correct X_1… X_6:

    ∏_k e^A^(k) X ∏_k'e^-A^(k') - e^A_c^(1) X e^-A_c^(1)≤  24· 2^0.5S(X)(24c'χ^ 56/225)^c+1X

Now, the A^(1)_i has support 4, but actually the first step of perturbation theory has support 3, just our definition overestimates it for the sake of consistent form over all k. So it's enough to group spins into pairs before applying the Trotter: A_μ,μ+1 = A^(1)_i + A_i+1^(1) and i=2μ-1. As a result we get an approximation e^-A_c^(1)_T by 4-local gates bounded using Lemma <ref> as:

    e^A^(1)_c X e^-A^(1)_c - e^A_c^(1)_T X e^-A_c^(1)_T ≤   O((4N+S(X)) A_μ,μ+1^2X/N) ≤ O((4N+S(X)) (c'χ)^2X/N)

and due to the grouping the condition of one circuit for the whole system is 4N≤ c. It can only be satisfied for c≥ 4, and N=floorc/4=1 always:

    ∏_k e^A^(k) X ∏_k'e^-A^(k') - e^A_c^(1)_T X e^-A_c^(1)_T  ≤  24· 2^0.5S(X)(24c'χ^ 56/225)^c+1X +O((4+S(X)) (c'χ)^2X)

Comparing the powers, we note that 56(c+1)/225>2 for c≥4, so the first term can be omitted for sufficiently small χ. We found that Trotter-approximating the first step of perturbation theory leads to an error

    ∏_k e^A^(k) X ∏_k'e^-A^(k') - e^A_c^(1)_T X e^-A_c^(1)_T  ≤ O((4+S(X)) (c'χ)^2X)

everywhere except for rare resonances. One can chose depth 2 and 4-local gates, or depth 3 and 3-local gates (the second application of Trotter bound leads to the same scaling of error). To go down to 2-local gates, one may apply Solovay-Kitaev to the 3-local gates, blowing up the depth by O(8^aln^b1/ϵ) = O(|lnχ|^b) factor. Note that even though for a uniform deep circuit the causal structure would make e^A_c^(1)_T X e^-A_c^(1)_T depend on O(|lnχ|^b) collar, the long patches of gates arising from Solovay-Kitaev are arranged in a checkerboard pattern, so the true causal cone remains within c.

To conclude, let us study the relationship between a state ∏_k e^A^(k) X|prod⟩ and e^A_c^(1)_T |prod⟩. Any observable X on S(X) sites is approximated by measuring it with e^A_T |0⟩ with error as in Eq. (<ref>). The latter measurement is represented by a system-size independent tensor network due to the cancellations. In particular, we can measure a matrix element of the density matrix over those S(X) sites. The difference in density matrices over S(X) sites is bounded element-wise:

    |(ρ -ρ_appr)_αβ| ≤   O((4+S(X)) (c'χ)^2X)

We can also show the operator norm bound:

    ρ -ρ_appr≤   O((4+S(X)) (c'χ)^2X)

here ρ = tr_S(X)∏_k e^A^(k) |prod⟩⟨ prod|∏_k'e^-A^(k') and ρ_appr = tr_S(X)e^A_c^(1)_T |prod⟩⟨ prod| e^-A_c^(1)_T. |∏⟩ is a product state in z-basis. To prove the norm bound, consider the eigenbasis of ρ -ρ_appr. We can make X to be measuring matrix elements in this basis. In particular, the maximum eigenvalue of ρ -ρ_appr will not be bigger than the difference |⟨ X ⟩ - ⟨ X ⟩_appr|≤   O((4+S(X)) (c'χ)^2X). Note that for the bound to be correct, the size for causal cone of X should be S(X) +2c and c≥ 4, and it does not improve the bound to have c>5.


To go to a better precision than (c'χ)^2, specifically χ^p, the depth required is exponential in p until we allow some threshold size S(X) +2c∼ p of gates in our circuit. Indeed, even for A = ∑_i A_i,i+1 above there isn't a subexponential way to achieve high precision using methods desribed so far. The threshold is reached when that ∑_j=1^c X_j needed to approximate to precision χ^p is contained entirely within one gate, so we can do ∏_k^k_max e^A_c^(k) in one step. This does not allow us to have a total circuit though - so methods described so far always require exponential in p size of the circuit if one wants to have one circuit working for operator in any place of the system that reduces to a collared one via cancellations. There might be a method of getting a total circuit of ∼ p local gates and subexponential depth via splitting e^A =(e^A_L⊗ e^A_R) e^G_LR where  L,R are half-chains one the left and right of a specific cut. The G_LR defined in this way does not commute with A_L, A_R, so the error above may be ∼A_R - divergent with the system size. In fact it's not! The G_LR has a collar bound/telescopic sum around the cut that can be found in <cit.>, but we do not pursue the construction of precise subexponential depth circuit here.












 §.§ Area law

We will prove the following:


 The state ∏_k=1 e^A^(k)|prod⟩ obeys the area law: for B - the half of the system, the entanglement entropy S(ρ_B) ≤ 650 c'χ. Here |prod⟩ is a product state in the z-basis.

We will also discuss the tail of the distribution of entanglement entropies in the presence of resonances.





We turn to the paper <cit.> which derives a bound on entanglement production S(t) -S(0) under the evolution e^iHt. In our case, we apply the bound for every step using iHt = A^k. We set t=1. 


We use the bound that is given on the last page of <cit.>. The interaction between two parts in H =-iA^k can be written as a Shmidt decomposition H_LR^(k) = ∑_n r_n J_L^n ⊗ J_R^n, where we have freedom to set ||J_L||=1, ||J_R||=1. The A^k_i are bounded by A_i^(k)≤ c' χ^ceilL_k-1 and there are S(A^(k))-1=floor(8/7L_k+1-1) of them across the cut, thus the interaction is bounded as H_LR^(k)≤ c'(S(A^(k))-1) χ^ceilL_k-1. We can decompose each A^k_i across the cut into 4^S(A^k_i,L) terms  — this is how much basis operators are there on one side of the cut. Each term in the decomposition is bounded as the whole operator that's decomposed A_i^(k):

    r_n≤ c' χ^ceilL_k-1

There are a total of ∑_S(A^k_i,L)=1^S(A^(k))/24^S(A^k_i,L)≤4/3 2^S(A^(k)) terms. The entangelement  produced over time t=1 will be 

    S(t) -S(0) ≤  c^* ∑_i r_i≤ c^*4/3 2^S(A^(k))  c' χ^ceilL_k-1

The coefficient c^*≈ 1.9 is universal and is found in <cit.>. The total entanglement is the sum of entanglement produced by each order k of the circuit. The total bound on entanglement is converging:

    δ S_tot≤∑_k  c^*4/3 c' 2^floor(8/7L_k+1)χ^ceilL_k-1≤
    ≤ c^*4/3 c'2^floor(8/7L_2)χ^ceilL_0∑_j=0^∞(2^floor(8/7L_3)-floor(8/7L_2)χ^ceilL_1-ceilL_0)^j ≤
    ≤ c^*4/3 c'2^7χ/1 -2^7χ≤  c^*4/3 c'2^8χ≤ 650 c'χ

where we replaced the sum of exp(L_k) =exp((15/8)^k)) by a geometric series with a factor such that every term upper-bounds the corresponding term in the original sum, and χ<1/2^8 has been used to get rid of the denominator. Thus the Imbrie states without resonances satisfy area law.

Every resonance is an extra local rotation over region of size n_r that does not have to be small. If it happens that our cut passes through the resonance, we need to add the entanglement created by this rotation to the total entanglement bound. A rotation over n_r sites can create Δ S(n_r) ≤ n_rln2 entanglement. Naive expectation is ≤ (n_r/2)ln2 which is true for a unitary acting on a whole system of size n_r. But in fact having a bigger system we can create two times more entanglement with the rotation on n_r sites!  Indeed, the bound is saturated by a swap gate,  that turns a state with n_r/2 Bell pairs on every side of the cut into a state where every Bell pair crosses the cut.
We believe that this is an overestimation and only two levels are mixed at a typical resonance, producing ∼ln2 entanglement.

After k_1 when the first resonance appears across the cut, there may be other, bigger resonances containing the first one. Since a cut is placed randomly in the system, denote the probability that the biggest resonance across the cut has size n_r by P(n_r). It is exactly the same quantity as the probability of a site to be a part of a connected cluster of resonances of size n_r. Imbrie shows that this probability decays faster than any polynomial. Thus the expectation value

    ∑_n_rn_rP(n_r) ≤const

Moreover, from investigating the first step of the construction we know that the first term has n_r=3 and P(n_r=3)=ϵ (with some coefficient depending on the probability distribution of disorder, in our case 1) and then it keeps decaying faster than any polynomial, so the leading contribution to the expectation value is the first term:

    ∑_n_rn_rP(n_r) ≤ 4ϵ

Here the first factor of 2 is the bound on entanglement instead of 3, as one of the sides has one site, and the second factor of 2 absorbs the later terms in the sum, which can be done for sufficiently small ϵ. Thus the total bound on the expectation value of the entanglement entropy is:

    S_tot≤   650 c'χ + 4ϵln2

We note that the first term γ^0.95 is smaller than the second term γ^0.05  — so the main contribution to average entanglement is from rare distributions of disorder when resonances cross the cut. We also note that in full construction by Imbrie, the perturbative parts of the circuit e^A^(k) get adjusted in the presence of resonances as the resonant spins are grouped into metaspins. The entanglement produced by resulting unitaries is still bounded, so we dropped this complication from our simplified circuit.

 Consider rare cuts that may have entanglement bigger than 2ln2, by passing through the big resonant regions. Assume for simplicity exponential  distribution of entanglement as a discrete random variable P(S=2xln2) = ϵ^x  for x the half of the size of the resonance. The biggest entanglement across the cut in a given system determined by:

    LP(S_max) = Lϵ^-S_max/2ln2≈ 1 
    
    S_max≈  2ln2ln L/(-lnϵ)

To derive maximal entanglement more carefully one needs to use the true statistics of resonances from Imbrie's paper, summarized in Postulate 4, that results in a superpolinomially decaying P(S) anyway. It only adjusts ln L → e^√(lnL) - a difference between the two will be almost impossible to detect.




 §.§ LRB




Consider 

    A(t) = e^iHtAe^-iHt =U^† e^iH_ztUAU^† e^-iH_ztU = ∑_r U^† e^iH_ztA_r e^-iH_ztU

where H is the Hamiltonian of our spin chain and U =∏_k=1 e^A^(k) is the unitary that diagonalizes it. According to Theorem <ref>, A_r supported on S(A) +2r sites and is bounded by 6· 2^0.5S(A)α^rA with α = 24c'χ^ 56/225. Now the terms in e^iH_zt that do not overlap with A_r commute through, so

    e^iH_ztA_r  e^-iH_zt = e^iH_z^r tA_r e^-iH_z^r t

where H_z^r = ∑_m H_z^r+m and the terms are bounded as 36α^m (using H=∑_i H_i, H_i≤ 2+γ≤3 and S(H_i) =2). Now we find m^* such that α^m^*t = 1. Define

    Q_r+m^* +p =e^i∑_m^m^*+pH_z^r+m tA_r e^-i∑_m^m^*+pH_z^r+m t- e^i∑_m^m^*+p-1H_z^r+m tA_r e^-i∑_m^m^*+p-1H_z^r+m t

and Q_r+m^* = e^i∑_m^m^*H_z^r+m tA_r e^-i∑_m^m^*H_z^r+m t. The higher-order Q's are bounded

    Q_r+m^* +p≤e^iH_z^r+m^* +p t - 1A_r≤Ac^*α^p

Where c^* is a new constant coefficient. We found that the telescopic sum for e^iH_z^r tA_r e^-iH_z^r t has the form

    e^iH_z^r tA_r e^-iH_z^r t = ∑_p Q_r+m^* +p,   Q_r+m^* +p≤A_rc^*α^p

After plugging it in the decomposition of A(t):

    A(t) =∑_r∑_p U^†  Q_r+m^* +p e^-iH_ztU

The last conjugation with the Imbrie circuit produces one more telescopic sum. Resumming all three telescopic sums, we get:

    A(t) =∑_q  Q_m^*',   Q_m^* +q≤A_r1/2 c^**α^q

A new constant and a factor 1/2 was introduced for convenience. Now a local operator B distance d_AB from A will only overlap with m^*+q ≥ d_AB terms in the telescopic sum. Resumming once more, we obtain:

    [A(t),B]≤ ||A|| ||B||c^**α^d_A,B-m^* = ||A|| ||B||c^**α^d_A,Bt


That concludes the proof.




§ INTERPRETATION OF IMBRIE'S PROOF



 §.§ Bound on generators


 Here we discuss the bound on the local generator A^(k)_i of the perturbative step e^∑_i A^(k)_i of Imbrie circuit. Here i is the spatial index labeling sites of the chain, and k is the order of the perturbation theory. Each A^(k)_i is acting on floor(15/7L_k) spins, where the radius of locality L_k = (17/8)^k is growing with k. We have presented a bound
 
    A^(k)_i≤ c'χ^ceilL_k-1

Imbrie works with a different decomposition: instead of A^(k) =∑_i A^(k)_i he uses A^(k)=∑_g A^(k)(g), where the sum iterator g will be explained below. On p. 27 of the arXiv version of the paper<cit.>, Imbrie presents a combinatorial argument that  ∑_g|i A^(k)(g) of terms contributing to A^(k)_i converges, and moreover the Bound (<ref>) holds, however the Bound is not explicitly presented.
 
Let's define the path formalism used by Imbrie in the following toy examples to facilitate the understanding of the argument on p. 27.



  
Toy example 1.
Consider an operator J defined as:

    J= ∑_g J(g)

The sum is over paths g = { i_1… i_|g|} in the space of spin flips. Each i_n is a site on a lattice, and the length of the string of those sites g is denoted |g|. In the ∑_g, it is implied that |g|>1. There's one condition on the string g: each new spin flip has to be within 1 lattice spacing from one of the old ones:

    ∀ n>1  ∃ n'<n : |i_n -i_n'|≤ 1

The operators J(g) are supported on the smallest consecutive region containing all spin flips in g plus a 1 spin collar on the left and on the right. Consider the following way to bound the terms in J:

    |J_σσ'(g)|≤δ^|g|/C(|g|)

where C(x) is a combinatorial factor to be determined and δ≪ 1. We'd like to see what C(x) do we need to be able to translate this bound to a more traditional norm bound. Consider a decomposition of J into "local" terms centered around sites of the lattice:

    J = ∑_i J_i,    J_i = ∑_g|i_1 =iJ(g)

where J_i contains only graphs starting at i_1. Let's show that a norm bound on J_i follows from a bound on J_σσ'(g):

    J_i≤∑_g|i_1 =iJ(g)

Note that any operator J(g) supported on ∪ g ∪ c where c is a 1-spin collar can be written in an operator basis with one basis vector per matrix element :

    J(g) = ∑_σσ'J_σσ'(g)|σ⟩⟨σ'|

and since J(g) acts as an identity outside ∪ g ∪ c, and doesn't flip c, the number of terms in the sum is 4^|∪ g|2^|c|≤ 2^2|g| +c:

    ∑_g|i_1 =iJ(g)≤∑_g|i_1 =i∑_σσ'J_σσ'(g)≤∑_g|i_1 =i2^2|g| +Cδ^|g|/C(|g|)

We can now split the sum into terms corresponding to each |g|:

    J_i≤∑_x =1^∞2^2x +Cδ^x/C(x)∑_g|i_1 =i,   |g|=x1

A simple upper bound on the number of elements in the sum ∑_g|i_1 =i,   |g|=x1 is if a flip i_n can be chosen from 2n-1 sites centered at i_1. For |g|=x flips, it gives (2x-1)!! = (2x-1)(2x-3)….

    J_i≤∑_x =1^∞2^2x +Cδ^x(2x-1)!!/C(x)

So a choice of C(x)=(2x-1)!! will cancel this factor. What remains is a geometric series which is bounded for 4δ<1/2:

    J_i≤ 2^C+1(4δ)

we have arrived at a desired bound.



  
Toy example 2. We will repeat the above construction with a minor modification: in the 

    J= ∑_g J(g)

the ∑_g is now taken over g such that |g|>L_1. The bound on J_σσ'(g) and C(x) =(2x-1)!! are the chosen the same. All the steps go through in exactly the same way, until we arrive to:

    J_i≤∑_x =L_1^∞2^2x +Cδ^x(2x-1)!!/C(x) = ∑_x =L_1^∞2^2x +Cδ^x≤ 2^C+1(4δ)^L_1

we see that any level is bounded by the smallest allowed flip.



  
Toy example 3. We will now introduce an extra layer of structure to match Imbrie's construction. Consider

    J= ∑_G J(G)

where G ={ g_1, g_2 … g_|G|} a string of spin flip paths g_n defined in examples before. The costraint on the paths in G is now:

    ∀ n>1  ∃ n'<n : d(g_n,g_n')≤ 1

where d(g_1,g_2) is the distance in lattice spacings between the two nearest spins in g_1,g_2. 

The support of J(G) is in the smallest consecutive region containing all spin flips of underlying paths of G plus some constant collar c. Surprisingly, it is enough to bound

    |J_σσ'(G)|≤1/C(|G|)∏_g∈ Gδ^|g|/C(|g|)

where C(x) = (2x-1)!! again. Indeed, for a "local" term

    J_i = ∑_G| i=i_1∈ g_1 ∈ GJ(G)

the norm bound

    J_i≤∑_G| i=i_1∈ g_1 ∈ GJ(G)≤∑_G| i=i_1∈ g_1 ∈ G4^∑ |g|2^c/C(|G|)∏_g∈ Gδ^|g|/C(|g|)

The sum over G can be split into sums over paths with fixed values of |G| and all |g_m|:

    J_i≤∑_|G|,|g_m|4^∑ |g_m|2^c/C(|G|)∏_mδ^|g_m|/C(|g_m|)∑_G| i=i_1∈ g_1 ∈ G,  |G|,|g_m| 1

so we need to do the counting again. We first count the small graphs as before, and then count how many ways are there to place |G| small graphs of sizes |g_m| so that each new one is a neighbor to one of the old ones. We get the following bound:

    ∑_G| i=i_1∈ g_1 ∈ G,  |G|,|g_m| 1 ≤
    ≤∏_m C(|g_m|) |g_1|(|g_1|+2|g_2|)(|g_1|+2|g_2| + 2|g_3|)… (2∑_m |g_m| -|g_1|)

or if we plug it in:

    J_i≤
    ≤∑_|G|,|g_m| 4^∑ |g_m|2^cδ^∑_m|g_m||g_1|(|g_1|+2|g_2|)(|g_1|+2|g_2| + 2|g_3|)… (2∑_m |g_m|-|g_1|)/C(|G|)

Or if we separate sums over |G|,∑_m |g_m|=X,|g_m|=x_m:

    J_i≤∑_|G|∑_X≥ |G| 4^X2^cδ^X∑_x_m| ∑_m x_m =Xx_1(x_1+2x_2)(x_1+2x_2 + 2x_3)… (2X-x_1)/C(|G|)


We see that the smallest power of δ is given by ∑_m |g_m| = |G|, in which case the C(|G|) exactly cancels its numerator:


    J_i≤∑_|G|  4^|G|2^cδ^|G| F(|G|,δ)

where

    F(|G|,δ)=1 +
     +∑_X> |G| 4^X-|G|δ^X-|G|∑_x_m| ∑_m x_m =Xx_1(x_1+2x_2)(x_1+2x_2 + 2x_3)… (2X-x_1)/C(|G|)

we want to show that  f(|G|,δ) can be bounded by a constant for sufficiently small δ. It is a series in terms of n=X- |G|. F = ∑_n=0^∞F_n and:

    F_n =(4δ)^n ∑_x_m| ∑_m x_m =|G|+nx_1(x_1+2x_2)(x_1+2x_2 + 2x_3)… (2(|G|+n)-x_1)/C(|G|)

There are always |G| terms x_m. One can already see that

    ∑_x_m| ∑_m x_m =|G|+nx_1(x_1+2x_2)(x_1+2x_2 + 2x_3)… (2(|G|+n)-x_1)/C(|G|)≤ (|G| +n)^2|G|

so we only get a polynomial prefactor in front of δ^n and the series converges. What we want to show is that  it converges to something that goes to ∞ only exponentially in |G|, so we can suppress it by (4δ)^|G|. Clearly the naive bound (|G| +n)^2|G| where dropped C(|G|) altogether gives a factorial-like scaling even for the first term in the sum that we know is 1. We know that number of terms in the sum for each n is in fact given by C_|G|+n^|G|< 2^|G|+n, so what's left is to bound a single term in the sum:

    x_1(x_1+2x_2)(x_1+2x_2 + 2x_3)… (2(|G|+n)-x_1)/C(|G|)≤
    ≤((2(|G|+n)-1/2|G|-1)^|G| = (1 +2n/2|G|-1)^|G|

these terms summed with 4δ^n can only lead to logarithmic corrections (1/(-lnδ))^|G| to the F(|G|,δ. Putting it all together:

    F_n ≤(4δ)^n 2^|G|+n(1 +2n/2|G|-1)^|G|

There is no factorial in |G|, so:

    F (|G|,δ) = ∑_n=0^∞F_n ≤ c'= O(1)

dropping the logarithmic corrections, and

    J_i≤∑_|G|c' 2^c (4δ)^|G|≤ c' 2^c+1 (4δ)

Each term in the sum is bounded by our choice of C(|G|). We can also consider long paths of |G|>L_2 like in the toy example 2 and get c' 2^c+1 (4δ)^L_2.



  
Actual construction


The structure of J_i that we described is present in the generators A_i of steps 1 and 2 of Imbrie's construction, with several differences. One difference is that in the original bound on A_σσ'(g) the γ/ϵ stands in place of δ that we used here. Another difference with J_i is that the paths for A are bounded in length: from below which leads to a higher power of γ/ϵ and from above by a growing scale L_k+1, which gives truly local generators, but does not in any way affect the derivation above. Also, the paths in g^(k) do not have to touch on a chain - instead, a collar is allowed, which is then lifted to fill in the gap with new graphs as described on pp. 45-46 of paper<cit.>. We discuss how the collar appears in the Subsection <ref>.

Finally, in Imbrie's paper C(x) (2x-1)!!∼ (x!)^1/2, instead it's an arbitrary power C(x) ∼ (x!)^2/9 (the graph factorial is defined on p.23 in the same way as our toy examples). The combinatorial factors are still controlled, but they require a more careful way to bound ∑_g and  extra tricks like resumming the long paths as a subroutine, as described on p. 27 of paper<cit.>. We note that in the earlier arXiv version of the paper, the n!^1/2 was used as a combinatorial factor.

Here we repeat the argument on p.27. The construction results in H
_pert=∑_g J(g). Then ⟨σ | H_pert |σ'⟩  =∑_g J_σσ'(g) where the terms in the sum have a bound proven in<cit.> for them. We'd like to prove that the sum is bounded as well. For that we do a transformation ∑_g |J_σσ'(g)| =∑_g |J_σσ'(g)| c(g) c^-1(g) ≤ max_gJ_σσ'(g)| c(g) ∑_g c^-1(g) where c(g) is some arbitrary function of the graph g. If c(g) is chosen in such a way that ∑ c(g)^-1≤ 1, then ∑_g |J_σσ'(g)| ≤ max_gJ_σσ'(g)| c(g). Imbrie chooses c(g) in such a way that max_gJ_σσ'(g)| c(g) ≤ (aχ) ^min|g| with some constant a, and proceeds to show that  ∑ c(g)^-1≤ 1 for that choice. Since we consider J_σσ'(g), it is nonzero only for g's that connect σ and σ', the min|g| is at least the length of the interval containing the difference σ-σ'. Similar process can be applied to A, where  min|g| is the minimal  remaining order ceilL_k-1.



 §.§ Bound on support



Imbrie's diagonalization proceeds in steps labeled by k=1,2,3…. At each step a unitary rotation O_re^A^(k) is applied to the Hamiltonian, where O_r is an operator acting spatially only near resonances, and as identity everywhere else. The generator of non-resonant part of the circuit is given as a sum of local terms A^(k) = ∑_g A(g), where g is the data in the format of the generalized path, containing, in particular, the locations of spin flips of A(g). The number of spin flips is upper bounded by the order of the generalized path |g|  — the power of the magnetic field involved in this term (some flips can act on the same location and "cancel" each other). At step k, that order is L_k-1≤|g| < L_k, where L_k= (15/8)^k (see, for example, p.24 of the new arXiv version). One may group terms into A^(k)_i = ∑_g:i∈suppA(g)A(g) where the first site on which A(g) acts nontrivially is i. We demonstated in Subsection <ref> how  A^(k)_i constructed this way can be bounded knowing the bound on A(g). Here we'll discuss the support of A(g) and A^(k)_i. 

The construction has a separate term A(g)^(k) for each matrix element A(g)^(k)∼ |σ⟩⟨σ'|, given by A(g)^(k) = J(g)^(k-1)/E_σ^(k-1) -E_σ'^(k-1). Here J^(k-1) are the terms appearing after k-1 previous rotations of H. From the point of view of support, J(g)^(k-1) combines support of the lower order paths that are part of g data, while E_σ^(k-1) is calculated using the diagonal part of the Hamiltonian truncated to terms J(g)_diag^(k-1) of order less than L_k-1 (see p. 48 of paper<cit.>. Such terms are diagonal, thus they have to wrap onto itself at least once to cancel all the flips. Using that, Imbrie proves a strong bound on extra support on each side c_k=1/2 suppE^(k-1)_σ - E^(k-1)_σ'≤15/14L_k-1. So the term we get in A(g) has the same positions of spin flips as J(g), and a collar around those spin flips of size 15/14L_k-1.



Let's now discuss the support of J^(k-1)(g) included in A^(k)(g). The orders of terms A^(k)(g) are in [L_k-1,L_k).

The order |g| limits the size the interval containing its spin flips. Thanks to bridging of the gap described on on pp. 45-46 of paper<cit.>, order directly translates into maximal size of the interval where flips are happening (it is exactly |g|<L_k for terms in A^(k), and |g|/2≤ L_k-1/2 for diagonal terms in J(g)_diag^(k-1) appearing in the denominator for A^(k)). So the total support is either |g| or |g|/2 plus two collars c_k. We illustrate how the collars appear by explicitly listing the properties of terms for the first few steps. 



 L_0=1     perturbation order [1,15/8] is removed at step 1 by A^1  ⇒     H^1 + J^1 
 
             J^1 starts with order 2 
  
 L_1=15/8     perturbation order [2,3] is removed at step 2 by A^2  ⇒     H^2 + J^2 

        H^2 includes order 2, support 3 terms     J^2 starts with order 2 
 
   L_2≈ 3.5     perturbation order [4,6] is removed at step 3 by A^3_c=2 ⇒     H^3 + J^3 

        H^3 includes order 6, support 7 terms     J^3 starts with order 7 
 
     L_3≈ 6.6     perturbation order [7,12] is removed at step 4 by A^4_c=5 ⇒     H^4 + J^4 

        H^4 includes order 12, support 16 terms     J^4 starts with order 13 
 
L_4≈ 12.4     c = 6+3 + 1+1 = 11         



We see that even though in H^3 the support 7 terms are present, somehow the collar for A^4 is chosen to be 5. In other words, the path in the collar at every step should be touching, otherwise the collars blow up. We assume that terms with bigger collar are not included in each A(g)^(k). Let us see what is the order of the extra terms that would appear if we tweak the procedure this way  — if it is the same order as J^(k), then we can just include them with other paths.


For illustration, let's use denominators within a shorter collar in A^4, and check that  e^A(H^3 + J^3) e^-A = H^4 + J^4 with the same condition on J^4 to include only |g|>13. Recall that we want to remove terms in J^3 of orders 7 to 12. Now in principle  H^3 contains order 6, support 7 terms, that look like 2-site collar around the graph of step 3. The collar itself is a step 2 graph of support 3, order 2, with collar of size 1 around it. It doesn't go further down. Now after conjugation, [A,H^3] should cancel terms in J^3 in such a way that only terms of order >12 remain. But instead of cancellation, it gives a (1- Δ E_full/Δ E_short)J^3_7..12 remainder Δ E_full -Δ E_short contains terms in H^3 whose center graph does not touch the spins flipped by A. In this case, it is only the center size 3, collar 2 (so support 7) term of order 6. 7+6 =13 which is enough to send it to J^4, and it's not even enough to make it connected - one needs to lift an order 2 term to make it continuous. Unfortunately this activity does not really fit into the structure that Imbrie proves bounds for, but it can be done at least in principle. Since Imbrie makes no mention of this, we leave the exact way it's done as an open question and proceed. 

Here's the table at k'th step:



     L_k-1 = (15/8)^k-1     pert. order [ceil(L_k-1),floor(L_k)] removed at step k by A^k_c_k ⇒     H^k + J^k 

        H^k includes order floorL_k, support floor(L_k/2)+ 2c_k terms     J^4 starts with order floorL_k+1 
 
L_k≈ (15/8)^k     c_k+1 = floor(L_k/2)+ c_k         



In other words c_k+1 = 1 +∑_k'=2^kfloor(L_k'/2) ≤∑_k'=1^kL_k'/2 which is the expression used by Imbrie. Let's evaluate it: 

    c_k+1 = 1/2∑_k'=1^k (15/8)^k' = 1/2(15/8)^k∑_k'=0^k-1 (8/15)^k'≤1/2L_k∑_k'=0^∞ (8/15)^k'  = 1/2L_k1/1 - 8/15 = 15/14 L_k

This is exactly as Eq 4.50 on p.46. So the support of a term in A^k is bounded by floorL_k + 2c_k ≤ L_k + 15/7 L_k-1 = 15/7L_k = 8/7 L_k+1. And the order is bounded by L_k-1 from below as seen in Subsection <ref> (even k=1 works out since ceilL_0=1). The final piece we need for the postulates in the main text is the support of resonances: after perturbative rotations A^k of step k, there come nonperturbative rotations that may flip spins within an interval 3L_k (p.15) and have a collar 15/14L_j = 15/14L_k-1 (p.25) on each side. So total support is 3L_k + 8/7 L_k = 41/7 L_k≤ 4.2 L_k.

