State Dependent Spread Complexity Dynamics in Many-Body Localization Transition
Maitri Ganguli Aneek Jana








[]maitrig@iisc.ac.in
Department of Physics, Indian Institute of Science, Bangalore




[]aneekjana@iisc.ac.in
Center for High Energy Physics, Indian Institute of Science, Bangalore





§ ABSTRACT 
We characterize the Many-Body Localization (MBL) phase transition using the dynamics of spread complexity and inverse participation ratio in the Krylov space starting from different initial states. Our analysis of the disordered Heisenberg spin-1/2 chain unravels that the ergodic-to-MBL transition can be determined from the transition of the pre-saturation peak in the thermofield double state (TFD) spread complexity. On the other hand, if an initially ordered state or a superposition of a small number of such states is chosen, then the saturation value of spread complexity and Krylov inverse participation ratio (KIPR) can distinguish the ergodic phase from the integrable phases, with no sharp difference between the integrable phases. Interestingly, the distinction between the disorder-free integrable and the MBL integrable phase is established by the spread complexity study of random states chosen from unitary and orthogonal Haar ensembles. We also study the complexity dynamics by coupling the system to a bath, which shows distinctive profiles in different phases. A stretched exponential decay of KIPR is observed when the MBL system is connected to the bath, with the decay starting at an earlier time for a greater value of environmental dephasing. Our work sheds light on the efficacy of Krylov space dynamics in understanding phase transitions in quantum many-body systems.


State Dependent Spread Complexity Dynamics in Many-Body Localization Transition
    Aneek Jana 0009-0001-1097-4250
    September 9, 2024
===============================================================================










  
*Introduction.—

In recent years, the study of quantum complexity quantified on the Krylov basis has gained significant interest for its usefulness in understanding the various
aspects of quantum many-body systems, quantum field theories, holographic models, quantum circuits etc <cit.>. The basic notion of complexity captures the difficulty of preparing a certain quantum state starting from a given initial state  <cit.>. In the context of quantum dynamics, Krylov Complexity measures the average position of a time-evolved state/operator in the Krylov basis formed by the action of the generator of time-evolution using the Lanczos algorithm or some modified version, such as the bi-Lanczos algorithm in the cases of non-unitary dynamics <cit.>. The Hamiltonian is rendered into a tridiagonal form in the Krylov space spanned by the Krylov basis vectors. Therefore, the complex quantum mechanical state/operator dynamics problem is effectively reduced to an equivalent single particle hopping problem in a semi-infinite lattice numbered by the index of Krylov basis vectors, where the hopping amplitudes are given by the Lanczos coefficients, which are the outputs of the Lanczos/bi-Lanczos algorithm.

While operator complexity in the Krylov space has been studied extensively, in a plethora of both closed and open quantum systems <cit.>, starting from the seminal paper by Parker et al.<cit.>, the study of the complexity of spread of states, also known as Krylov spread complexity, is relatively new <cit.> and its relevance has been investigated in various contexts including integrability to chaos transitions <cit.> and 𝒫𝒯-symmetric non-Hermitian Hamiltonians <cit.>, etc very recently. This work aims to contribute to this endeavor by studying state-dependent spread complexity dynamics in the systems that exhibit Many-Body Localization Transition (MBLT) [Krylov operator complexity in MBL is studied in <cit.>].

The discovery of many-body localization (MBL), where strong disorder and interaction lead to emergent integrability, is a great example of the violation of Eigenstate Thermalization Hypothesis (ETH) <cit.> beyond integrable systems. In the interacting systems, the presence of disorder <cit.> or quasiperiodicity <cit.> generically can give rise to MBL (which is a generalization of Anderson localization <cit.> to interacting systems); however, in thermodynamic limit and greater than one dimension its stability is a subject of active debate <cit.>. Recent experimental findings have provided direct evidence of this breakdown of ergodicity in interacting many-body systems. Specifically, these studies have observed such behavior in various systems involving ultracold atomic fermions <cit.>, a chain of trapped ions <cit.>, and also in superconducting circuits <cit.>. In these systems, strong disorder caused them to become localized, thereby preventing them from reaching thermal equilibrium as expected in the absence of such disorder.

Though for both MBL systems and disorder-free integrable systems, the existence of an extensive number of conserved quantities is observed, which gives rise to the absence of level repulsion (hence Poisson level spacing distribution in the energy spectrum <cit.>), but from the point of view of the entanglement entropy there is a difference. The MBL systems possess eigenstates showing area-law entanglement entropy, whereas generally, with few exceptions, disorder-free integrable systems show volume-law entangled eigenstates. So, from this perspective, identifying the MBL system from thermal eigenstates is easier than identifying the integrable systems. The localization behavior of MBL prevents the system from exploring all possible states, effectively breaking ergodicity—the principle that a system will eventually explore all accessible microstates if given enough time. We are interested in exploring this behavior of MBL systems through the spread complexity and differentiate the MBL emergent integrability from disorder-free integrable systems in this perspective.

 Recent studies have focused on the spread complexity of the Thermo-Field Double (TFD) state, which is a canonical purification of the Gibbs density matrix, to distinguish chaotic and integrable phases <cit.>, as well as to observe the integrable-to-chaotic transition by treating the peak in the spread complexity as an order parameter<cit.>. Our study verifies that the peak in the TFD spread complexity indeed shows a transition in the ergodic-to-MBL crossover. 
However, in the present work, our scope is more focused on the ergodic-to-MBL transition in a more elaborate way by studying the initial state dependence of the spread complexity in various physically relevant scenarios. We comment on how the distinguishability among the phases can be seen in the dynamics of initial states chosen from random Haar ensembles, making a clear distinction between the strong-disorder MBL integrable phase and the integrable phase that exists in no or weak disorder in finite-size systems [we should mention that in the strict thermodynamic limit, even an infinitesimal disorder can break integrability, see <cit.>]. 

Finally, coupling the system to a bath weakly, we have demonstrated that the MBL system shows stretched exponential decay of the Krylov space localization measure (similar to the decay profile of an initial density pattern <cit.>).

Apart from the Krylov spread complexity, we have also looked at the Krylov Inverse Participation Ratio (KIPR) to understand the dynamics of various states in the Krylov space clearly. The KIPR gives a good dynamic measure of how a state is localized on the Krylov basis. 
All these complexity measures are ultimately dependent upon the wave-function coefficients of the state in the Krylov basis; nevertheless, they allow us to shed light on different aspects of the dynamics in the Krylov space.






  
*The Model and the Method.—

In this work, we have considered the paradigmatic model of Many-body Localization <cit.><cit.>, the spin-1/2 Heisenberg model with random-field disorder,


    H = 1/2∑_i (σ_i^xσ_i+1^x+σ_i^yσ_i+1^y + σ_i^zσ_i+1^z) + ∑_i h_i σ_i^z


here the random-fields h_i are sampled from a uniform distribution [-W,W]. It is known that by increasing the disorder strength, the system goes through an MBL transition around W ≈ 3.5 <cit.><cit.>. This is an instance of a chaotic-to-integrable transition. Another disorder-free integrable phase exists (at W=0) and is sustained at a very small value of W for finite size systems <cit.>. We have tried to characterize both of  these integrable and the ergodic phases through the analysis of Krylov complexity, starting from various initial states. We work in the zero magnetization sector (∑_i σ^z_i=0) for definiteness and use the periodic boundary condition.

The first complexity measure we use is the Krylov Spread Complexity (KSC). Starting from an initial state |ψ_i⟩, we calculate the orthonormal Krylov basis vectors {|K_n⟩} (using Lanczos algorithm for Gram-Schmidt orthogonalization on the set {|ψ_i⟩,H|ψ_i⟩,H^2|ψ_i⟩,…}) and then expand the time-evolved state in this basis,

    |ψ(t)⟩ = e^-i H t|ψ_i⟩ = ∑_n ϕ_n (t) |K_n⟩


The ϕ_n(t)'s can be caculated numerically or using the Lanczos coefficients {a_n} and {b_n} (see Supplemental Material <ref>) in solving the following recursive differential equation,

    i d/dtϕ_n(t)= a_n ϕ_n(t)+b_n ϕ_n-1(t)+b_n+1ϕ_n+1(t),

with the boundary condition given by ϕ_n(0)=δ_n,0.

If we denote the probability of being in the n-th basis vector by p_n then,

    p_n(t) = | ϕ_n (t) |^2,   ∑_n p_n(t) = 1


The Krylov spread complexity is given by the average position in the Krylov basis,


    𝒞_𝒦 (t) = ∑_n n p_n (t)


While the above complexity measures average position, we need something that measures the typical number of basis elements needed for describing the time-evolved state by an entropic notion. For this purpose we use the Krylov entropy,

    𝒮_𝒦(t) = - ∑_n p_n(t) log p_n (t)


and associated Krylov Entropic Complexity (KEC),

    𝒞_𝒮 (t) = e^𝒮_𝒦 (t)


As a measure of localization in the Krylov space, one can define the Krylov Inverse Participation Ratio (KIPR),

    ℐ_𝒦(t) = ∑_n p_n^2(t) ≤ 1


Where larger values of ℐ_𝒦 will imply localization in Krylov space. ℐ_𝒦 also satisfies the lower bound ℐ_𝒦≥ 1/d_𝒦≥ 1/d, where d is the dimension of Hilbert space under consideration and d_𝒦 is the dimension of Krylov space with d_𝒦≤ d.


 We also consider the effect of dissipation in this system due to its coupling to the environment via a thermal bath or some measurement apparatus. To do so, we have used an effective non-unitary evolution of the state of the system, specifically, the evolution of a single quantum trajectory, which corresponds to the no-click/no-jump limit in a suitable post-selection procedure.
 
For simplicity, we consider only two jump operators that couple to the system with a coupling strength or dephasing α. The jump operators are chosen to commute with the total magnetization operator to keep the state in the same magnetization sector <cit.>,

    L_1 = σ_0^xσ_1^x + σ_0^yσ_1^y
       L_2 = σ_L-2^xσ_L-1^x + σ_L-2^yσ_L-1^y

The non-Hermitian Hamiltonian is,

    H' =  H - iα(L_1^† L_1+L_2^† L_2)


To elevate the notion of spread complexity for non-Hermitian evolution, one needs to use the bi-Lanczos algorithm, which reduces to the usual Lanczos algorithm in Hermitian limit <cit.>. Here we have two sets of Krylov basis vectors which are bi-orthogonal, {|P_n⟩} and {|Q_n⟩}, and we can expand the time-evolved state in these two sets as following,

    |ψ(t)⟩ = ∑_n ϕ^q_n(t) |P_n⟩ = ∑_n ϕ^p_n(t) |Q_n⟩

For this scenario, one can modify the notion of probability by introducing additional necessary normalization,

    p_n (t) = |(ϕ^p_n(t))^*ϕ^q_n(t)|/∑_m |(ϕ^p_m(t))^*ϕ^q_m(t)|

Once the probability of being in the n-th Krylov basis is defined, the definitions of the various complexity measures are kept unchanged.

In the following sections, we emphasize the complexity dynamics for various initial states and the effect of dissipation in the presence of different disorder strengths W. Further, we study the distinctions between the phases that exist for different ranges of W.





  
*Thermofield Dynamics.— To understand the complexity dynamics, at first, we have chosen the Thermofield Double (TFD) state, which is an entangled state in the product Hilbert space of the two copies of the same system. If the system has a spectrum {E_n}, and corresponding eigenvectors |Ψ_n⟩,  then the TFD state at inverse temperature β is defined by,

    |TFD(β)⟩ = 1/√(Z_β)∑_n e^-β E_n/2|Ψ_n⟩_L ⊗|Ψ_n⟩_R

where Z_β is the thermal partition function and L and R denotes left and right copies in the product Hilbert space, respectively. For this state's time evolution, one can consider only the time evolution of the left copy by ℋ, and the right copy does not evolve. So effectively, we do the dynamics of the following Gibbs state,

    |ψ_β⟩=1/√(Z_β)∑_n e^-β E_n/2|Ψ_n⟩


By choosing β=0, we have done the exact time evolution of this state and the corresponding complexity 𝒞_𝒦, with different choices of W, which gives a characteristic peak in the 𝒞_𝒦 evolution for the chaotic regime, and the peak disappears in the MBL regime.

The plot of peak height (here we define peak height by (max(𝒞_𝒦(t)/d)-0.5), which we take as an order parameter) against disorder strength for various system sizes (see Fig. <ref>) helps us to determine the range of disorder where the transition occurs, which agrees with earlier estimates from level-statistics transition <cit.>.

From Fig. <ref> one important point to note is that the early growth of complexity is controlled by W monotonically (this can be related to the fact that the minimum energy difference in the spectrum increases with W, <cit.> shows that the time-scale in spread complexity growth depends on the minimum energy difference in the spectrum), even though for larger W, the complexity fails to reach the peak.



 The other two complexity measures shed light on a different aspect of the dynamics in the Krylov space. We observe that (see Fig.<ref>) both the integrable phases, disorder-free integrable and MBL integrable, delocalizes faster in the Krylov space at early times than the deep ergodic phase (W≈ 1). However, the KIPR and the entropic complexity, which measure localization in Krylov space, are the same in the ergodic and MBL phases at late times. In contrast, in the disorder-free integrable phase, the final state is significantly localized in comparison. Understanding the behavior of early time localization in Krylov space of the TFD state can be an important direction for further research.

We have successfully probed the chaotic-to-integrable transition in the context of MBL transitions using the complexity dynamics of the TFD state. However, one should remember that the dynamics of TFD states do not use the information of the associated eigenvectors; it solely depends upon the distribution of eigenvalues. So, studying the dynamics of the TFD state alone will leave us with an incomplete picture of the actual MBL transition. This motivated us to consider the complexity dynamics of other states, which we discuss in the following sections.








  
*Relaxation Dynamics.—



It is known that MBL systems break ergodicity and can retain some local information about the initial state  <cit.> . Since our complexity measures essentially give us the information on how distant the time-evolved state is from the initial state and its spread in the Krylov basis, it does make sense to use the spread complexity of states as a probe to measure this memory-retaining property of the MBL phase and its initial state dependence [In <cit.>, the authors studied memory-retaining property of MBL phase for different ordered initial states using the statistics of Lanczos coefficients.].

For the model under consideration, we prepare the initial state in the computational basis, which has a Néel-like order [The authors in <cit.> considered the domain-wall-like state as the initial state and found different complexity behavior in ergodic and MBL phases. Our results are more general and explain their observations.],

    |ψ_i⟩ = |1010⋯10⟩

where 1 (0) at ith position denotes the eigenstate of local σ_i^z operator with eigenvalue +1 (-1). The results discussed below hold whenever the number of computational basis elements on which the initial state has support, is much less than the total Hilbert space dimension. We show the complexity dynamics for the initial state in Eq.(<ref>) in Fig. <ref>.

Fig. <ref> shows that the late-time averaged value of all three complexity measures point towards the fact that the MBL system retains significant memory of the initial state (Robustness of this behavior of retaining the initial memory will be probed by coupling the system to the environment, in a later section). On the other hand, the ergodic phase makes the state more complex and delocalized than its integrable neighbors, as clear from its higher complexity and lesser KIPR.

The above observations can be explained by the existence of pseudo-spin-like quasi-local integrals of motions (LIOMs) <cit.> or local-bits (ℓ-bits) τ_i^z, which have finite overlap with the local spin operators σ_i^z in the presence of strong disorder. Since LIOMs are conserved quantities, information encoded in their initial values remains intact unless the system is coupled to a bath. Therefore, the computational basis elements, eigenstates of the local spin operators σ_i^z, show less complex and more localized dynamics in the MBL phase. This is the same reason why, in the MBL phase, the operator Krylov complexity of local σ_i^z operators show more localized behavior (in Krylov space) <cit.>.






 Now, if we consider extensive superpositions of computational basis elements as our initial state, then the characterization of different phases through complexity becomes more involved. To capture the typical behavior of complexity dynamics in the different phases, we need to be more careful in choosing initial states. To probe typical behavior, we choose Haar random states from two different ensembles, unitary and orthogonal, as our initial states, which we discuss next. 





  
*Complexity of Typical States.— We have been choosing various initial states from the angle of different physical motivations. We have observed that some states, for example, the TFD state and initially ordered states, carry the direct signatures of integrability, be it the disorder-free integrable or MBL integrable phases. Despite these successes, one should try to understand the complexity dynamics of typical states in different phases to see whether the distinctions among the phases via complexity dynamics are generically present.

To be completely generic at first, we choose states distributed randomly in the N-dimensional complex projective space ℂℙ^N according to the Haar measure (where (N+1) is the dimension of the Hilbert space). Such states can be sampled by the action of random (N+1)× (N+1) unitary matrices on some chosen arbitrary state. We then compute the complexity dynamics of such random Haar states evolving under specified Hamiltonian. 
We find that the MBL phase can be distinguished from the ergodic phase from the absence of a peak in the complexity profile and absence of a dip in the KIPR profile. On the other hand, the disorder-free integrable phase can be distinguished from its lower complexity and higher KIPR, see Fig. <ref>.





Another choice can be to sample random states uniformly from the real projective space ℝℙ^N, which can be physically important while being completely random. Such states can be sampled by acting random orthogonal matrices on some arbitrarily chosen state. It is again observed that a peak in the complexity can distinguish the ergodic phase from the integrable ones. Even though, at late times, all of them have similar saturation values, unlike the ℂℙ^N. But from Fig. <ref>, it is clear that the KIPR value still can set apart the disorder-free integrable phase from the MBL phase.

From these observations, we infer that the distinctions among the chaotic and the integrable phases are not special to states like TFD; rather, they occur in more generic states also. Therefore, it is worth understanding the physical origin of such distinctive complexity behavior of different kinds of random states under the evolution of different Hamiltonians.











  
*Dissipative Dynamics.—
MBL is recognized as a robust dynamical phase of matter, when strictly decoupled from the environment, does not thermalize <cit.>. However, coupling to the environment necessarily leads to a delocalization transition, with a time scale governed by the coupling parameter value. The goal is to understand this transition from a Krylov complexity perspective and point out the robustness of the MBL integrability.

  We model the coupling to a bath or a measurement apparatus by an effective non-unitary evolution by the Hamiltonian in Eq.(<ref>), which describes a specific quantum trajectory. Under the assumption of weak dephasing (α<<1), this effective description is justified, and as we will demonstrate, it well captures the essential physics even within this minimal setup. To be specific, we start with the Néel state, whose complexity dynamics in a closed system can distinguish the ergodic from the MBL, to see how its evolution is affected by the environmental coupling.
  
   Spread complexity and KIPR for different disorder strengths W are affected by non-zero α, which is shown in Fig. <ref>. It is found that coupling to the environment certainly causes delocalization in the Krylov space and an increase of spread complexity. We have observed that, by increasing disorder W, the early-time growth rate of spread complexity increases first and decreases by further increasing of W. In particular, in the large W phase, the early-time growth rate is comparably less [However, in the late-time spread complexity behavior, we have found that in the large W phase, the complexity is more than the chaotic and disorder-free integrable phase. To understand the precise reason, further investigation is needed. We would like to address this point in a future work.]. Therefore, the MBL phase is less susceptible to environmental coupling from the Krylov complexity perspective, at least at early times.
  
  Now, focusing on the delocalization transition in MBL systems for different values of the environmental coupling α, we have observed that the Krylov IPR has a stretched exponential decay profile (see Fig. <ref>), which is similar to the decay profile of initially set particle density imbalance found earlier in <cit.> and decay of third Renyi negativity as found in <cit.>[In this context, one should note that the idea of KIPR is more robust in the angle of understanding the localization property for the system. We do not need to use Néel like order to observe this kind of characteristic decay of KIPR in open quantum systems, in fact, any computational basis state would capture this feature.].
  However, the complexity increases always for non-zero α (see Fig. <ref>) with a possible saturation above 0.5 at very late times. This indicates that the time-evolved state really goes far from the initial state in the Krylov space at late times. These observations show that the quasi-local integrals of motion (LIOMs) are destroyed when the MBL system is coupled to a bath, which is efficiently captured by Krylov space dynamics.

  
  
  

  

  




  
*Discussion and Future Directions.—
In our work, we have studied the state-dependent spread complexity dynamics in the disordered spin-1/2 Heisenberg chain to understand whether complexity dynamics in the Krylov space can make distinctions among the disorder-free integrable, ergodic, and emergent integrable MBL phases. 


    
  * Through our analysis we show that the pre-saturation peak height in TFD state complexity, which we use as an order parameter, can significantly capture the ergodic to MBL transition.
    
  * If we start with an initially ordered state like Néel state we can also infer different phases from the saturation values of Krylov Complexity and Krylov Inverse Participation Ratio (KIPR).
    
  * Our work has established that not only special states like TFD or initially ordered states but also randomly chosen typical states carry important information about the phases in their complexity dynamics.
    
  * However for the open system if we focus on the early-time growth of the complexity of the initial Néel state, in the MBL phase it is slower than the other two phases. Interestingly, for the MBL phase coupled to a bath, we observe the dissipative delocalization in the Krylov space with a stretched exponenetial decay profile.

   
So, our analysis highlights the effectiveness of Krylov space methods in understanding the different phases for both closed and open quantum dynamics.

One important future direction is to do the complexity dynamics in the quasiperiodic systems that show MBL transition, such as the interacting Aubry-Andre model <cit.>, to understand if there is any qualitative difference in the complexity dynamics of disordered systems and the quasiperiodic potential systems. Investigating the role of interaction in the complexity dynamics is also interesting. In the spin model, this amounts to tuning a parameter that multiplies the term σ_i^zσ_i+1^z. One should also try to understand if time-reversal symmetry (TRS) has any role in this context. For that, a term that breaks time-reversal symmetry <cit.> is to be added to the Hamiltonian. Then, the complexity dynamics should be studied for this Hamiltonian that breaks TRS.

Finally, we leave it for future work to analyze the spread complexity dynamics in the MBL phenomenological model in terms of the LIOMs <cit.>. That can potentially provide a more model-independent way of characterizing MBL transition through the spread complexity dynamics.








M.G. would like to thank Sumilan Banerjee and Subroto Mukerjee for useful discussions. The authors also thank Aranya Bhattacharya for comments on the draft and suggestions. M.G. is supported by the Integrated PhD fellowship of Indian Institute of Science, Bengaluru, and A.J. is supported by INSPIRE fellowship by DST, Govt. of India.








  
Note added.– After the completion of this work, we became aware that the authors <cit.> are investigating complexity in random unitary circuits, which also have signatures of MBL like behavior from a Krylov complexity perspective.

apsrev4-1











§ SUPPLEMENTAL MATERIAL




§ BI-LANCZOS ALGORITHM AND COMPLEXITY


The bi-Lanczos algorithm is suitable for tri-diagonalizing a non-Hermitian matrix (H^†≠ H). This can be applied for non-Hermitian state evolutions (in non-unitary dynamics, no-click limit in MIPT or no-jump limit in open quantum system evolution) or operator evolution (in the context of open quantum system dynamics and general measurement settings). In the context of state evolution, the Hamiltonian H will be non-Hermitian, and in the context of operator evolution, the Lindbladian ℒ_o (ℒ_o^†≠ℒ_o) will be non-Hermitian. Once the inner product is specified, the following bi-Lanczos algorithm can be used for both non-Hermitian state and non-Hermitian operator evolution. For Hermitian generators, this reduces to the original Lanczos algorithm. This algorithm is based on <cit.>. 


Suppose we trying to tri-diagonalize an operator M which is non-Hermitian, that is M^†≠ M. In this case, we have to construct two sets of Krylov basis vectors {|P_n⟩} and {|Q_n⟩}, which are bi-orthogonal, 

    ⟨ Q_m | P_n ⟩ = δ_mn


So these two sets are orthogonal with respect to each other, but inside each set, the vectors are not orthogonal; by that we mean ⟨ Q_m | Q_n ⟩≠δ_mn and ⟨ P_m | P_n ⟩≠δ_mn.


As we will see, we shall generate three sets of Lanczos coefficients, the main diagonal {a_n}_n≥0, the upper diagonal {b_n}_n≥1 and the lower diagonal {c_n}_n≥1. We shall have b_n = c_n for the Hermitian limit, but they would be different for non-Hermitian.



  
Constructing the Krylov basis vectors. Start with |P_0⟩ = |Q_0⟩ = |ψ(0)⟩ or |𝒪_0⟩ for state and operator complexity respectively. And define the initial Lanczos coefficients values a_0 = ⟨Q_0| M |P_0⟩,  b_0 = 0,  c_0 = 0. The following algorithm is similar to the usual Lanczos algorithm where the |P_n⟩'s are constructed using M and |Q_n⟩'s are constructed using M^† but in a correlated way.


For n=0, 

    |A_1⟩ = M |P_0⟩ - a_0 |P_0⟩
       |B_1⟩ = M^†|Q_0⟩ - a_0^* |Q_0⟩
        w_1 = ⟨A_1|B_1⟩,   c_1 = √(|w_1|),   b_1 = w_1^*/c_1
       |P_1⟩ = |A_1⟩/c_1,  |Q_1⟩ = |B_1⟩/b_1^*
        a_1 = ⟨Q_1| M |P_1⟩


For n≥ 1, 

    |A_n+1⟩ = M |P_n⟩ - a_n |P_n⟩ - b_n |P_n-1⟩
       |B_n+1⟩ = M^†|Q_n⟩ - a_n^* |Q_n⟩ - c_n^* |Q_n-1⟩
        w_n+1 = ⟨A_n+1|B_n+1⟩,   c_n+1 = √(|w_n+1|),   b_n+1 = w_n+1^*/c_n+1
       |P_n+1⟩ = |A_n+1⟩/c_n+1,  |Q_n+1⟩ = |B_n+1⟩/b_n+1^*
        a_n+1 = ⟨Q_n+1| M |P_n+1⟩


It can be shown that the two sets of Krylov basis vectors as formed by above algorithm are indeed bi-orthogonal. In practical purposes one has to stop once c_n+1 is less than some cut-off and one has to full-orthogonalize |A_n+1⟩ in the |Q_m⟩ basis and full-orthogonalize |B_n+1⟩ in the |P_m⟩ basis for m≤ n.


Observe that,


    M |P_n⟩   = a_n |P_n⟩ + b_n |P_n-1⟩ + c_n+1|P_n+1⟩
    
            M^†|Q_n⟩   = a_n^* |Q_n⟩ + c_n^* |Q_n-1⟩ + b_n+1^* |Q_n+1⟩


It shows that M is rendered to a tridiagonal form in the |P_n⟩⟨Q_m| basis.




  
Complexities.

Now expand |ψ(t)⟩ or |𝒪 (t)⟩ in the both sets of Krylov basis vectors (let's call them P-basis and Q-basis vectors).


    |ψ(t)⟩ = ∑_n ϕ_n^q (t) |P_n⟩ = ∑_n ϕ_n^p (t) |Q_n⟩


where,


    ϕ_n^q (t) = ⟨Q_n|ψ(t)⟩,   ϕ_n^p (t) = ⟨P_n|ψ(t)⟩


For non-Hermitian evolution we define a probability P(t) by,


    P(t) = ∑_n |(ϕ_n^p (t))^* ϕ_n^q (t)|


Let's define the Krylov complexity as average position in the Krylov basis dictated by the wavefunction in this basis,


    𝒞_𝒦(t) = ∑_n n |(ϕ_n^p (t))^* ϕ_n^q (t)|/∑_m |(ϕ_m^p (t))^* ϕ_m^q (t)| = ∑_n n p_n


where we have the diagonal probabilities,

    p_n = |(ϕ_n^p (t))^* ϕ_n^q (t)|/∑_m |(ϕ_m^p (t))^* ϕ_m^q (t)|


The Krylov entropy can be defined as the following Shannon entropy of the diagonal probabilities,


    𝒮_𝒦(t) = - ∑_n p_n log p_n


Then we define entropic complexity as exponential of the Krylov entropy,


    𝒞_𝒮(t) = e^S(t)


We also define the Krylov Inverse Participation Ratio (KIPR),


    ℐ_𝒦(t) = ∑_n p_n^2


The KIPR signifies the extent of localization of the wave function in the Krylov basis, and the Krylov entropic complexity signifies the effective number of Krylov basis vectors on which the wave function has typical support.




§ LEVEL STATISTICS AND STATISTICS OF LANCZOS COEFFICIENTS

The transition of level statistics is a useful probe for chaotic-to-integrable transition. For chaotic or ergodic systems, we have the Wigner-Dyson level statistics, and for integrable systems, we have the Poisson level statistics. However, during chaotic to integrable transition, the level statistics change smoothly from Wigner-Dyson to Semi-Poissonian to completely Poisson. A flow can model this whole transition in the space of different random matrix ensembles, specifically Gaussian β-ensembles, by changing the value of β (note that β=1 corresponding to Gaussian Orthogonal Ensemble (GOE) and β=0 correspond to Poisson distribution) <cit.>, see Fig. <ref>.



The variance of Lanczos coefficients also captures chaotic and integrable behavior. In particular, integrable systems will have a higher variance of Lanczos coefficients, and chaotic systems will have comparatively less variance of Lanczos coefficients, see Fig. <ref>, Fig. <ref>. This is observed in the case of operator growth also. In the context of the spread complexity of TFD state in random matrix theory (RMT) ensembles, the variance of Lanczos coefficients has anti-correlation with the average level-spacing ratio <cit.>. This holds for our model also; see Fig. <ref> for the initial TFD state and Fig. <ref> for the initial Néel state.












§ EIGENSTATE IPR AND SPECTRAL FORM FACTOR


One of the important characteristics of the many-body localized phase is the localization of many-body eigenstates in real space due to the presence of disorder. This can be quantified by evaluating the inverse participation ratio (IPR) of many-body eigenstates |Ψ_n⟩ in the computational basis elements |i⟩, IPR = ∑_i |⟨Ψ_n|i⟩|^4. A higher value of IPR will capture the localization of the eigenstates. In Fig. <ref> we have plotted the IPR of many-body eigenstates for different disorder strengths. It is observed that the disorder-free integrable phase has the most delocalized eigenstates; on the contrary, the MBL phase has the most localized eigenstates. Therefore, we observe that the localization of eigenstates increases monotonically with the disorder strength. 

While the eigenstate IPR sheds light on the localization property of the MBL phase, it does not capture the emergent integrability aspect. Emergent integrability can be understood from the level-spacing ratio as well as spectral form factor (SFF). Spectral form factor is defined by,

    SFF(t) = 1/d^2∑_n,m=1^de^i(E_n-E_m)t


For chaotic eigenvalue distribution, the SFF shows a characteristic dip-ramp-plateau behavior as a function of time t. However, this feature is absent in integrable models. From Fig. <ref>, we observe that the ergodic phase SFF has indeed a dip-ramp-plateau behavior, which is absent in both the disorder-free integrable phase and the MBL phase.









§ QUENCH DYNAMICS

we consider some possible quench scenarios to understand state-dependent spread complexity in different phases. To do so at first we assume the disorder-free phase to be at its ground state, then suddenly we add disorder in the system with varying strengths. Moreover, we track the complexity dynamics of the state under the evolution of the disordered Hamiltonian. Our observation indicates that the complexity and delocalization both increase with increasing strength of disorder, and the normalized spread complexity seems to saturate to around 0.16 towards high disorder (top-panel in Fig. <ref>). In contrast, if one starts with an infinite temperature TFD state corresponding to the disorder-free model and evolves with disordered Hamiltonian, then in the high disorder limit, the normalized complexity seems to saturate to around 0.36 (bottom-panel in Fig. <ref>).



Further, to understand the reversed quench scenario, we take the initial state, which is the ground state of a phase inside the MBL regime, say W=6.5, and evolve it with a Hamiltonian with lesser disorder strength. At W=6.5, it will have zero complexity since the complexity evolution of eigenstates is trivial. However, for lesser disorders, the state will evolve non-trivially, and we observe that the complexity saturation increases with the decrement of W. In our case for system size L = 12 as we go away from W=6.5, complexity increases with a maximum saturation around 0.4, but as we enter the disorder-free phase, the complexity saturation decreases to around 0.29 for W=0 (middle-panel in Fig. <ref>).

Taking into account the observations and previous results on thermofield dynamics, relaxation dynamics and quench dynamics, it gives an indication that the disorder-free phase generally exhibits lower complexity compared to the other two phases. The results for typical state spread complexity in the main text also support this observation.








