Operator growth and Krylov Complexity in Bose-Hubbard Model

We study Krylov complexity of a one-dimensional Bosonic system, the celebrated Bose-Hubbard Model. The Bose-Hubbard Hamiltonian consists of interacting bosons on a lattice, describing ultra-cold atoms. Apart from showing superfluid-Mott insulator phase transition, the model also exhibits both chaotic and integrable (mixed) dynamics depending on the value of the interaction parameter. We focus on the three-site Bose Hubbard Model (with different particle numbers), which is known to be highly mixed. We use the Lanczos algorithm to find the Lanczos coefficients and the Krylov basis. The orthonormal Krylov basis captures the operator growth for a system with a given Hamiltonian. However, the Lanczos algorithm needs to be modified for our case due to the instabilities instilled by the piling up of computational errors. Next, we compute the Krylov complexity and its early and late-time behaviour. Our results capture the chaotic and integrable nature of the system. Our paper takes the first step to use the Lanczos algorithm non-perturbatively for a discrete quartic bosonic Hamiltonian without depending on the auto-correlation method.


Introduction
Quantum mechanics is a staple ingredient in the playground of every theoretical physicist and is also responsible for most modern developments in natural sciences.A basic idea during our introductory lessons in quantum mechanics is the evolution of a simple operator O with time t for a system governed by a Hamiltonian H.In the Heisenberg picture, the simple operator O at t = 0 grows in the Hilbert Space into a complex operator O(t) through e iHt Oe iHt .This growth quantifies the spread of local information supported by a single or small number of sites to the entire system.The nature of the operator growth depends on the system's Hamiltonian.An operator, for example, will grow into an integrable and non-integrable system in distinct ways.Using the Baker-Campbell-Hausdorff formula, the operator growth is expressed in terms of complicated nested commutators [H, [H, We soon realize that the evaluation becomes rapidly difficult to perform with each level of the nested commutator.Nevertheless, this simple notion of operator growth plays a massive role in diverse areas of physics.The 'dissipative' behaviour of the operator growth is recently found to be connected to the thermalization of a given system [1][2][3][4][5].
Another common way of describing the operator growth in the space-time picture is by the out-of-time-ordered (OTOC) commutator [1,2,[6][7][8].The classical chaos makes the OTOC grow exponentially, also known as the Lyapunov exponent.The interesting conjecture is that the Lyapunov exponent is bounded from above [9] for a few strongly coupled large-N models (the Sachdev-Ye-Kitaev (SYK) model, for example, having a dual in the gravity side) [10][11][12].However, the bound on the Lyapunov exponent is not seen universally, especially for generic, non-integrable Hamiltonians outside of the semi-classical, large-N limit [3,7,13,14].
A general theory of operator growth under generic, non-integrable Hamiltonian dynamics is recently given by Parker et al [15].The authors argued that, with each increasing level of the nested commutators in the operator growth, the growing operators increase exponentially in time.The exponential size of the computation enables the system to be essentially treated as a thermodynamic entity, thus, emerging a statistical picture.This gave rise to the idea of governing the operator growth with some universal properties, even for non-integrable systems in arbitrary dimensions.The hypothesis for universal operator growth in [15] is based on the notion of a recursive technique, the Lanczos algorithm.Starting with a test operator |O(t = 0)), the Lanczos algorithm iteratively generates a set of orthonormal basis operators.This basis of orthonormal operators is called the Krylov basis |O n ), which can have fewer elements than the dimension of the Hilbert space, depending on the dynamics and the choice of the initial operator.Along with the basis operators |O n ) (where n = 0, 1, • • • , denotes the dimension of Krylov space), a set of coefficients {b n } is generated, known as Lanczos coefficients.The time evolved operator |O(t)) is expanded in terms of the Krylov basis |O n ) and time-dependent coefficients ϕ n (t), which obeys discrete Schrödinger equations [15].The operator dynamics is encoded in these time-dependent coefficients, which can also be thought of as wavefunctions distributed over the Krylov basis.The operator growth can be understood as the delocalization of an initial operator which we choose to be |O 0 ) over the Kyrlov basis with time, and the hopping amplitudes are given by the Lanczos coefficients.
The concept of various notions of complexity follows naturally from the above observation.For example, Parker et.al [15] have introduced the notion of Krylov complexity, (or K-complexity) as a measure of quantum chaos in the thermodynamic limit, that depends on the local Hamiltonian, a specific choice of an inner product, and the choice of the reference operator.K-Complexity of a generic operator in a finite-dimensional system size exhibits a specific profile with time which was first shown in [16] and then in [17].With an aim to find the chaotic nature of the system, we perform a numerical analysis to find the profile of K-Complexity at all time scales, including the scrambling phase and the late time phase, and compare our results to those of [16,17].
The Lanczos coefficients b n , in a quantum chaotic system, are hypothesized to grow fast.The maximum possible growth varies linearly with n [15], modified by a logarithmic correction in local (1 + 1)-dimensional lattice systems.This is in contrast, to the integrable models, [16] where the Lanczos coefficients grow as b n ∼ αn δ , where 0 < δ < 1.This behaviour of the Lanczos coefficients b n , implies that the K-complexity grows exponentially in time in chaotic systems and powerlike in an integrable system.In recent times this has motivated lots of studies.Interested readers are referred to some of these references, which are by no means exhaustive .
Motivated by all these, we have initiated a study of operator growth for the Bose-Hubbard model in this paper.It has found important application in the field of ultracold atom physics [61].This model exhibits Superfluidity-Mott insulator phase transition [62], which has been observed experimentally with ultracold atoms in optical lattice [63].More importantly, this model provides us with a rich parameter space, tuning by which it can behave as chaotic or non-chaotic (integrable) systems [64][65][66].Although for the system with fewer degrees of freedom, this kind of mixed behaviour can be expected [67][68][69], but for a generic quantum many-body system mixed systems are less common [70].So this model provides us with an ideal playground to study operator growth as well as the associated Krylov complexity for both the chaotic and integrable regime and make a comparison.In [70], the authors have considered a classical limit by parametrizing the interaction by β = U N , where U is the onsite interaction and the particle number is N , and then showed that in the limits β → 0 and β → ∞, the three-site Bose-Hubbard model possess integrable nature due to the vanishing of the only positive Lyapunov exponent λ max .Also, this model presents us with an excellent opportunity to study the operator growth and the Krylov complexity for a quantum-many body system consisting of bosonic degrees of freedom system containing quartic interaction terms using the Lanczos algorithm non-perturbatively.To the best of our knowledge, in the past, only limited studies of this kind have been made for non-quadratic systems 1 .
The paper is organized as follows.In Sec. 2, we introduce the necessary details of the celebrated Bose-Hubbard model, used in the context of our current analysis.In Sec. 4, we dive into the details of the Krylov basis and spanning of the Krylov space, which is significantly different from the Hilbert space.Then, in Sec. 5, we revisit the Lanczos algorithm.Even after being a powerful technique, Lanczos method suffers terribly from instabilities and errors rounding up in the process.This lead to this Sec.5, where we discuss the resolution to apply the modified version of the Lanczos method.Our results of the operator growth for the Bose-Hubbard model is presented in Sec.6 and that of the Krylov complexity is presented in Sec. 7. Finally we end with a summary of our results and future directions in Sec. 8. We also have incorporated Appendices (A), (B) and (C) explaining more details of our numerical analysis. 1Non-Gaussian RMT are studied in [57,71], but they are 0 + 1D systems.For SYK see [15,18].However, they are fermionic systems.For studies in JT gravity see gravity see [30] and for holographic CFT see [25].Also, in [72] a study of Krylov complexity has been made for λϕ 4 upto first order in λ.Unlike our case, authors have used the auto-correlation method [15,25] to compute Lanczos coefficients and they work in the continuous limit.Also for a study of the Krylov complexity for Calabi-Yau quantum mechanics please refer to [36].But this is a single particle quantum mechanical system.

Bose-Hubbard model
In this section, we start by summarising the necessary details of the one-dimensional Bose-Hubbard model.The Bose-Hubbard model consists of interacting bosons placed on a lattice.The model rose to prominence to describe the dynamics of ultracold bosonic atoms in an optical lattice, with the system parameters being completely controllable by laser light [61].We restrict the model to open-boundary chains of length M with a vanishing chemical potential.The Hamiltonian is given by where < i, j >≡< j, i > denotes the summation over adjacent sites (j = i ± 1), a † i and a i are the bosonic creation and annihilation operators on the site i (with the site i running from 1 to the site number).Also, ni = a † i a i is the number operator that counts the number of particles on that site i.Here, M denotes the site number, and N is the particle number.In the first term, the summation runs up to M − 1 when we consider the fixed boundary condition, such that particles are to be arranged in a linear chain with non-interacting walls at the end.Another possible boundary condition is to put the particles in a circle with the M + 1-th site being identified as the first site.In this periodic boundary condition, the summation over the first term runs up to M .J is the tunnelling coefficient which signifies the energy gained due to the particles hopping from one site to its neighbouring site.In this paper, we will only use open boundary condition.The first term proportional to J is the kinetic term of the Hamiltonian.On the other hand, the parameter U characterizes the two-particle on-site repulsive interaction strength.A point of importance is that the parameters J and U can be adjusted by various means.The model with a finite value of chemical potential shows two different phases depending on the ratio of the two parameters Λ ≡ U J measured in terms of the energies.For Λ << 1, the phase is identified by weak coupling and strong hopping, giving rise to the super-fluid phase.Additionally, Λ >> 1 with strong coupling and weak hopping characterizes the Mott-insulator phase [73].
The Hamiltonian has a U (1) symmetry associated with the conservation of a total number of particles.The symmetry group of the open boundary Hamiltonian can therefore be written as a U (1).For closed boundary conditions, the symmetry group of the Hamiltonian becomes U (1) ⊗ D M , where the symmetry D M is the combination of translation and reflection symmetries.The dimension of the Hilbert space H with the total number of atom N and number of sites M is found to be 2) The dimension D grows with the number of particles N as D ∼ N M +1 for a fixed M .A classical limit exists in the Bose-Hubbard model when the particle number N → ∞ with the site number M is kept fixed.The classical limit is similar to the discreet non-linear Schrodinger equation.The interaction is usually parametrized in the classical limit by β = U N .The model in the β = 0 and β = ∞ limits are integrable, thus, analytically solvable.
For the site number M = 2, the model is solvable for all values of β.Otherwise, for M ≥ 3 and finite β, the Bose-Hubbard model is non-integrable [70]2 .However, for M = 3, the model shows neither strongly chaotic nor integrable behaviour.It is said to exhibit highly mixed behaviour.Our findings in Sec 6, also confirm the said behaviour.In the following sections, we analyze the Bose-Hubbard model for site number M = 3 with different values of the particle number N and Λ ≡ U J .It is possible to generalize our results for higher particle and site numbers.However, the computations grow excessively expensive and are thus limited by the computational power.We will study the time evolution of an operator with a fixed Hamiltonian Eq.(2.1) in Krylov space.We use the Lanczos algorithm to find the appropriate Krylov basis and the corresponding Lanczos coefficients.To the best of our understanding, this is the first time the Lanczos algorithm (appropriately modified, discussed in later sections) is being used for a bosonic lattice system.

Revisiting the Krylov Space
In this section, we review the details of the Krylov Basis and its relation to the Krylov space.The Krylov basis is an ordered, orthonormal set of base kets that contains information about the time evolution of an operator.We begin by choosing an appropriate basis associated with the Hilbert space H with dimension dim(H) = D.The associated Hilbert space of linear operators L(H) ≡ Ĥ that acts on H, is of dimension D 2 .Next, we choose a system of N particles in M sites governed by the Hamiltonian (H) defined in Eq.(2.1) and an operator O (H, O ∈ Ĥ) [17].A natural choice of basis for the Bose-Hubbard Hamiltonian is the base ket with the occupation numbers in each site, i.e.
Here, n i ≥ 0 and N i=1 n i = N .One can construct all the basis vectors with appropriate N , and M following the rules mentioned in [74].
Next, we choose an initial operator O ∈ Ĥ and evolve it in time.The operator O gradually explores a subspace of the space of operators of the given system during its evolution.This subspace is called the Krylov space.The orthonormal Krylov basis K and the structure of the Krylov space can be analyzed using the powerful algorithm known as the Lanczos algorithm [15].The Lanczos algorithm also captures the rate at which a typical operator span the Krylov space.The Krylov basis K depends on the choice of the initial operator, and the dynamics of the system may not expand the entire Hilbert space H.We use the smooth-ket notation |O) to denote the element of the operator space corresponding to an operator O throughout the paper.In the following paragraph, we explain the notions of the Krylov Basis and Krylov space in the context of an operator growth O(t).
In the Heisenberg picture, the time evolution of |O) is given by Here, L is the Liouvillian operator defined by L|O) = [H, O] with H being the system Hamiltonian Eq.(2.1).Thus, using the Baker-Campbell-Hausdorff formula in Eq.(3.2), the Krylov space H O associated with the operator O, can be defined as the linear span of all nested commutators of the Hamiltonian with the operator: The dimension of the Krylov space dim(H O ) ≡ K, is less than the dimension of the operator space Ĥ, (dim( Ĥ) = D 2 ).For most systems, even though this set mentioned in Eq.( 3.3) has infinite elements, only a few of them are linearly independent.One can do a spectral decomposition [17] of O in the energy basis in order to obtain a precise result : where {E a , |E a ⟩} D a=1 are the eigenvalues and eigenstates (respectively) of H.It can be shown following [17] that the elements of the sum in (3.4) are the eigenstates of the Liouvillian L : with |ω ab ) ≡ |E a ⟩ ⟨E b |, and (phases) ω ab = E a − E b .In general, the Krylov dimension K is equal to the number of nonvanishing projections of O over the eigenspaces of the Liouvillian [17].Therefore, the upper bound of the dimension K is given by : This consists only of the unavoidable degeneracies, i.e, the zero eigenvalue phase ω aa = 0 and has the degeneracy of D. The upper bound K = D 2 − D + 1 is only saturated for cases with no degeneracies other than the diagonal null phases ω aa .The dimension of the Krylov space decreases with the increase of degeneracies in the Liouvillian spectrum.
Ideally, a typical Liouvillian operator in a chaotic system has no degeneracies other than the null phases and is known to saturate [17] the upper bound of Eq.(3.6) while for integrable systems, the Krylov dimension K is substantially smaller than the bound.

Details of the Lanczos method
After reviewing the necessary details of the Krylov space, we are now ready to construct it for our system.In this section, we review the Lanczos algorithm [75,76] to construct the Krylov Basis and, subsequently, the Krylov Space before discussing its shortcomings.Given a certain inner product in the operator space, one can construct an orthonormal basis for Krylov space associated with an operator O, with a fixed Hamiltonian H.This operation is done by the Lanczos algorithm, which is another version of the Gram-Schmidt orthogonalization procedure.In this paper, we make use of the infinite-temperature inner product, namely the Frobenius inner product3 : where ||.|| is the induced norm.The Lanczos algorithm generalizes to the following systematic steps [17,75,76]: |A n ) ̸ = 0 for all n < K as the orthogonal basis K has a rank K or K independent directions.At the next iterative level, , which is already a complete orthonormal basis of H O and must therefore vanish.This explains the reason for b K = 0 in Step 5 and the termination of the iteration.Therefore, we conclude that the Lanczos algorithm must terminate once all the independent directions in H O are exhausted.
One important aspect of writing the operator space H O in the Krylov basis K is that the Liouvillian L in the Krylov basis simplifies to a tridiagonal matrix (O m |L|O n ) = L mn .The entries of this tridiagonal matrix are given by the Lanczos coefficients.Alternatively, The Lanczos sequence {b n } behaves in a particular manner for chaotic and integrable systems.In both cases {b n } initially grows and then decreases to terminate for finite-dimensional systems.The nature of the {b n } sequence is fluctuating.Therefore, in [22] the authors have introduced the variance of Lanczos coefficients as a possible relation between the chaoticity of the system and it is given by, where i runs from 1 to K − 1 4 , and ⟨• • • ⟩ represents the mean value.The variance σ 2 depends on the nature of the b n distribution.Before we end this section, we should keep in mind that, even after being a very powerful algorithm, the general Lanczos algorithm discussed above fails to construct an orthonormal basis often due to numerical instabilities.These numerical instabilities can occur due to the accumulation of errors at each level of the iteration.These instabilities can be removed using the reorthogonalisation technique discussed in the next section.

Reliability of Lanczos algorithm?
One of the most significant consequences of using the Lanczos algorithm is that it is supposed to ensure the orthogonality of the Krylov basis obtained through each iteration in n.The algorithm is designed to stop when the Lanczos co-efficient b n hits zero.However, the Lanczos algorithm suffers major instabilities in numerical calculations [17,[77][78][79].While computing the Krylov elements |A n ) in step 3 of the Lanczos algorithm described in the Sec. 4, we need the information of its previous two elements.Hence, the errors arising from the finite precision keep piling up, resulting in significant overlaps between the Krylov elements in each iteration.The more iterations, the faster the orthogonality is lost between the Krylov basis.
This results in the Lanczos coefficients b n growing initially but never reaching zero.The b n 's start to oscillate drastically around some average value after the initial growth.This unfaithful behaviour of the Lanczos coefficients was also reported in complex SYK [17], and our analysis of the Bose-Hubbard model also went through similar errors.
Nevertheless, there exist several ways to deal with the numerical errors piling up in the Lanczos algorithm, namely Full-orthogonalization (FO) and Partial re-orthogonalization [17,77,80,81].We found a useful resolution for our case, by implementing the Fullorthogonalization procedure. 5The key to implementing the FO algorithm is to execute the Gram-Schmidt at each iterative step n twice and also use the information about all the Krylov elements (not only the previous two, as described in the previous section) to confirm the orthogonality between the Krylov basis (upto the error coming due to the machine precision).It is beneficial to explicitly mention the FO algorithm we have used successfully for the Bose-Hubbard Model.We start by computing the zeroth-operator |O 0 ) from the initial operator |O) using Eq.(4.1) and step 2 of the Lanczos sequence.
3 Perform re-orthogonalization using all of the previous Krylov elements by

Numerical results
Now armed with the above discussion, we present our numerical analysis of operator growth for the Bose-Hubbard model, the Hamiltonian of which is given by Eq.(2.1).In particular, we study operator growth in Krylov space and the associated Lanczos coefficients respectively for a three-site and three-particle model, (M, N) = (3, 3), (3,5), and for a (3, 7) model with different values of Λ ≡ U J .We restrict ourselves to two ranges of values of Λ: Λ → 0, and Λ = finite where we expect that the integrable and chaotic behaviour will appear following our previous discussion in Sec. 2.
First, we choose an appropriate hermitian operator for numerical computations.The operator that we choose is given by: This choice of the operator has both the number operator at j th site nj and the imaginary hopping term.One can, in general, also choose other operators, such as the number operator at a particular site ni = a † i a i or an imaginary hopping term i(a † j a j+1 − a † j+1 a j ).We found that for all the choices mentioned above of the initial operator, the behaviour of the Lanczos coefficients b n , and Krylov complexity C K typically remains the same 6 and the conclusions drawn in this paper remain same.One can further refer to Appendix (A) where we have shown the numerical results for a different choice of initial operator (A.1).The operator growth in the Krylov basis for the initial operator mentioned in Eq.(6.1) using the Lanczos algorithm is shown for the first two steps in the Appendix (C).Below, we discuss the numerical analysis done to find out the full Lanczos sequence for different system configurations.We restrict our analysis to the a three-site Bose-Hubbard model with particle numbers N = 3, 5, 7. We state our numerical results for the three-site Bose-Hubbard model with an initial operator mentioned in Eq.(C.2), and for different values of interaction strength parameter Λ as following: • For Λ = 0, the Lanczos iteration stops at the fifth iterative order with b 5 = 0. So, the dimension K = 4 of Krylov basis is much smaller than the dimension of the Hilbert space of operators D 2 = 100 and hence, confirms the bound mentioned in Eq. (3.6).For smaller values of Λ ∼ 10 −7 − 10 −3 , the Lanczos iteration stops at the fifth iterative order 7 with b 5 ∼ 10 −6 .This kind of behavior is expected as the Bose-Hubbard model is learned to be integrable [70] for Λ → 0 as Λ = 0 or U = 0 is just a Gaussian model and hence this is consistent.8 • Similarly, for large Λ ∼ 10 15 , the Lanczos iteration9 stops at a very small value of K, which is in agreement with the integrable nature of the Bose-Hubbard model at Λ → ∞ [70] as in this limit, the only possible Lyaponuv exponent is found to be exactly zero.
• The Lanczos iteration saturates the upper bound mentioned in Eq.(3.6) for finite values of Λ.The plots of b n with n show a phase of initial growth up to n ∼ ln D. The initial phase of growth is then followed by a regime of steady decrease towards zero with a generally constant negative slope of the order ∼ −10 −2 ∼ K −1 .This is shown in the Fig. ( 1a), (1b) and (1c) for different values of (M, N ).
These observations suggest that the operator growth via the nature of the Lanczos coefficients can distinguish the system as being chaotic or non-chaotic (integrable), providing further evidence to the earlier claims [15,22,43].We summarize our conclusions regarding the nature of the Lanczos coefficients in the Table 1.

Variance of Lanczos coefficients
Before we end the section, we numerically study the variance10 of b n due to their fluctuating behaviour, for the integrable and chaotic regimes11 as a function of the parameter Λ.  (3,5).It can be seen from the plot that as Λ becomes finite starting from Λ = 0, the variance decreases and becomes close to zero for Λ = 1.This behaviour of the variance points towards the transition from integrable to chaotic.
• In Fig. 2, we have shown the nature of the variance σ 2 as a function of the interaction parameter Λ.We notice that as the value of the parameter Λ increases from 0 to 1, the variance σ 2 becomes closer to zero.The above behavior of the variance is typically true for chaotic systems [22,50] and hence we conclude that the system is integrable for Λ = 0 and becomes more and more chaotic as Λ approaches the value 1.
• Also, as we increase the system size, the Lanczos sequence smoothens out for the chaotic system, bringing the variance closer to zero and causing a more smooth transition from the integrable regime, Λ = 0 to the chaotic one Λ > 0.
• One can notice a sharp fall-off in σ 2 values as the value of the Λ increases from zero.It marks the transistion of the system from a integrable to a chaotic one as was expected from our previous analysis of the Lanczos coefficients.
• Also, one can find that a similar sharp peak can be obtained at Λ ∼ O( 102 ) which marks the transition from chaotic to an integrable system. 12 Also, in Fig. 2 one can notice "kinks" near Λ ∼ 0.1, Λ ∼ 0.15 etc.These kinks are due to the unsmooth behavior of Lanczos coefficients for small system size.As we increase the system size and b n s become smooth, these kinks disappear.One can immediately get a hint of this from the N = 5 plot where some of the kinks disappear.
We have listed our findings from the numerical analysis of b n pointwise above where we have analysed three-site Bose-Hubbard model for different system configurations.Note that the disordered average of b n has been considered in [17] for the SYK 4 model and in [27] for a many-body localization (MBL) systems, e.g., one-dimensional (1D) spin chain with L sites and open boundary conditions, where in both cases values for the interaction strengths are drawn randomly from a certain random distribution and then {b n } are averaged over all the random realizations.The interaction parameter Λ in the Bose-Hubbard model does not follow from a random distribution.We do our analysis for particular values of Λ to characterize the chaotic and non-chaotic (integrable) behaviour.

Krylov Complexity
Now we investigate the Krylov complexity, often termed as K-complexity that was introduced by Parker et.al [15] as a useful diagnostic of quantum chaos in thermodynamic limit which measures an operator's complexification over the orthonormal Krylov basis.[15].Also motivated by our study of the variance of Lanczos coefficients in the previous section, we will mainly focus on the finite Λ regime where the model is expected to show chaotic behaviour.It is obtained by first constructing the Krylov basis using the Lanczos algorithm given by (5) and then using the Lanczos coefficients b n to reduce the analysis of the time-evolution of an operator O into the solution of a differential recurrence equation [15].The definition of Krylov complexity is outlined below.The time-evolved operator can be expanded in the Krylov basis as: where ϕ n (t) can be thought of as single-particle wavefunctions distributed over the Krylov basis, which via the Heisenberg equation dO dt = i[H, O] satisfies the following differential recurrence equation, To solve this equation we need an initial conditions which is: For a normalised operator O(0) = O 0 , and ϕ −1 (t) ≡ 0 to make Eq.(7.2) consistent with Eq. ( 7.1) we need to impose this initial condition.Furthermore, the boundary conditions b 0 = b K = 0 are fixed to ensure the finiteness of this operator-evolution problem.One can also infer from the equation of motion mentioned Eq. ( 7.2) of the wavefunctions ϕ n (t) as a singleparticle hopping problem on a semi-infinite chain called the Krylov chain, with the hopping amplitudes given by the Lanczos coefficients b n and ϕ n (t) are the wave-packets traveling on the chain [15].Since the initial operator is normalized at the first step of the Lanczos algorithm mentioned in Sec. ( 5), the notion of unitarity suggests that the wavefunction ϕ n (t) is normalized at all times: Then Krylov complexity or K-complexity is defined as the time-dependent average position of the distribution ϕ n (t) over the ordered Krylov basis [15]: We discuss the behavior of the K-complexity emerging from the numerical analysis in the following points: • We observe an exponential growth in K-complexity C K (t) at very early times until t s = log (log D).This is known as the scrambling time.The value of scrambling time t s changes as the system configuration changes which is shown in Fig. 3. • After the scrambling time t s , the K-complexity C K (t) features a linear increase in time until times of order t ∼ D and finally at very late times of order t ∼ D 2 it saturates (with fluctuations around that saturation value as shown in Fig. 4) typically at the values 13 of order K 2 .This phenomenon can be understood physically as a uniform distribution of operator in Krylov basis.The behavior of the K-complexity at all time scales is shown in Fig. 4 for different values interaction parameter Λ.
The above results are summarised in Table 2.The characteristics of K-complexity described above clearly agree with the results of [16] where the profile of K-complexity for generic operators in a finite-size system was first shown 14 .One can refer to Appendix.(B)where numerical fit is performed to analyse the behaviour of K-complexity at different time scales.

Discussion
In this article, we considered the Bose-Hubbard model as a quantum many-body model that is neither integrable nor strongly chaotic.It also has a quartic term in the Hamiltonian, providing a platform to study Krylov complexity beyond the quadratic discrete bosonic quantum many-body system.The ultra-cold atom Hamiltonian is directly related to the Bose-Hubbard Hamiltonian.We investigate a considerably more basic form of the model with a chemical potential µ = 0 and, thus, does not exhibit any phase transition.We start by analyzing the operator growth in Krylov basis in a three-site (M = 3) Bose-Hubbard model for different parameter Λ values.Furthermore, we perform a numerical analysis to explore operator growth in the Bose-Hubbard model.We study cases with different particle sizes, N .Although we have shown the results for a specific initial operator, overall features remain the same even if we consider different initial operators.An important observation of our work is that the standard method of performing the Lanczos algorithm fails for this model due to the piling of errors at each level of the iteration.We remove this difficulty by performing the Full orthogonalization (FO) method to find the appropriate Krylov basis and the Lanczos coefficients.
One of our crucial findings was that the dimension of the Krylov basis (K) is much smaller than the dimension of the Hilbert space for the operators when Λ → 0 or Λ → ∞.This confirms the fact that the three-site Bose-Hubbard model becomes integrable at the asymptotic limits of Λ.The Lanczos coefficient b n displays ramp and plateau behaviour for finite values of Λ.This indicates that for finite values of Λ, the three-site Bose-Hubbard model becomes chaotic, and the dimension of the Krylov basis saturates the upper bound as mentioned in Eq. (3.6).Also, in addition to the study of the behaviour of Lanczos coefficients b n , we also study the variance σ 2 of the Lanczos sequence as a function of the interaction parameter Λ which arises naturally as {b n } are fluctuating in nature.We observed that the variance follows the typical nature of being close to zero [50] as we approach from an integrable system to a chaotic one.It will be interesting to figure out in future how we can extract the information about spectral statistics, which encodes the information of the quantum chaos from the statistics of b n perhaps following the study of [50,57].Furthermore, we investigate the Krylov complexity for finite Λ.We have investigated the behaviour of K-complexity at three different time scales.We discovered that the K-complexity exhibits exponential increase at initial time scales, typically of O(log(ln D)), followed by a linear growth over time.K-complexity saturates in time at very late times, of the order of the dimension of the operator Hilbert space.This demonstrates that the three-site Bose-Hubbard model becomes chaotic for finite Λ.The emergence of these three distinct time scales matches with the general claim made in the literature, e.g.[16].
There are various future directions that will be interesting to pursue.One can study the Liouvillian spectrum and operator matrix elements in energy eigenbasis to find universal chaotic behaviour of the system from ETH approach [85].Throughout our studies, we have set the chemical potential to zero.However, the system can undergo quantum phase transitions in the presence of non-vanishing chemical potential.It will be interesting to perform our analysis for this case and whether we can comment on the phase transition using Krylov complexity.This will provide a possibility of comparing with the results coming from the circuit complexity analysis for the Bose-Hubbard model as well as connecting with results from holography [86][87][88].Also, it will be interesting to extend our computation for finite temperature case and compare with the results of [45,72,89] whether we get a staggeredlike behaviour for Lanczos coefficients.Recently emergence is scar-like states which evade thermalization has been demonstrated in Bose-Hubbard type models [90][91][92].It will be interesting to study the (spread) complexity of such states along the line of [51].Furthermore, we like the continuous limit of this model study the Krylov complexity and hope to connect with the recent studies of it for field theories [25,72,89].We hope to address some of these issues in the near future and report on them.Last but not least, as mentioned earlier, the Bose-Hubbard model emerges from various physical scenarios involving ultra-cold atom gas in periodic potential (e.g.dimerized magnets [93], interaction-induced tunnelling [94]) and can potentially be used as a quantum-simulator [95], we hope that our study of operator growth and the associated complexity for this system might shed important light in this context and initiate further studies.

A An alternate choice of initial operator
In this section, we provide the numerical results obtained from studying the Lanczos coefficients and Krylov complexity for a choice of the initial operator that is different from the initial operator considered in the main text.We work in M = 3 and N = 5 system.The appropriate hermitian operator that we choose for numerical computations is given by the number operator at the first and second sites: The numerical results obtained are shown in the Fig ( 5).One can take other choices of different hermitian operators that do not commute with the Bose-Hubbard Hamiltonian (2.1).
For instance, the choice of the initial operator where the number operator at all sites is taken commute with the Hamiltonian (2.1) and thus makes the Lanczos algorithm trivial.

B Numerical fits of K-complexity at different time scales
Below we discuss the numerical fits done to obtain the different behaviours of C k (t) at different time scales.We have shown a case where Λ = 0.1, N = 5, and M = 3.We found that the exponent α = 0.49635 controlling the exponential behaviour of K-complexity at early times is approximately equal to the slope of b n , β = 0.53538 in the linear regime as was expected.where the subscript n represents the site number.We start with the three-particle and threesite Bose-Hubbard model where Λ = 2.5, and the initial operator is given by where n1 = a † 1 a 1 and n2 = a † 2 a 2 are the number operators on the 1st and 2nd site respectively.
• O 0 : The operator O ini is not normalised.We normalise it using the inner product defined in 4.1 to get our starting operator O 0 as: The first Lanczos coefficient is: The matrix representation of the operators A is written in the D = 10 dimensional Hilbert space.As a result, at this level of iteration, the orthogonalized Krylov operator is: Note that, the operator O 0 is a double-site operator.It grows to O 1 , which is a mix of double-site and quadruple-site operators.At the next step, the operator O 1 will grow to a six-site operator contribution to give O 2 .However, the operators cannot grow more than a six-site contribution because of the restriction on the total number of particles N = 3.Thus, the operator remains as a six-site operator until it exhausts all the linearly independent contributions.The rest of the orthonormal operators of the Krylov basis can be generated using the same approach.In this case, the algorithm terminates at the n = 92 level of iteration giving b 92 = 0. Hence the Krylov dimension is K = 91 which exactly saturates the bound mentioned in 3.6.For a larger value of N , the operator growth does not stop at a six-site operator.Hence, the analysis becomes more and more tedious and can be performed by a numerical analysis owing to the results in the main text.

4
with b n = ||A n || , 5 If b n = 0 stop; otherwise set |O n ) = 1 bn |A n ) and go to step 3.The above algorithm constructs an orthonormal basis for H O , |O n ) K−1 n=0 known as the Krylov basis and a set of coefficients b n K−1 n=0 known as the Lanczos coefficients.Thus, starting from |O 0 ), each iteration of the Lanczos algorithm produces an element |A n ) orthogonal to all previous Krylov basis elements |O m ), with m < n.
key difference from the usual Lanczos algorithm) , 4 Perform step 3 again (another key difference from the usual Lanczos algorithm) , 5 Now, define the Lanczos coefficient b n = (A n |A n ) , 6 Next, set Krylov basis |O n ) = 1 bn |A n ) and return to step 2 , 7 Stop when b n = 0 .

Plot of bn vs n for different values of Λ with N = 3 , and M = 3 ( a )Plot of bn vs n for different values of Λ with N = 5 , and M = 3 ( b )Plot of bn vs n for different values of Λ with N = 7 , and M = 3 (Plot of bn vs n for Λ = 2 . 5 , M = 3 , and N = 10 ( d )
Lanczos sequence for M=3, N=3.Lanczos sequence for M=3, N=5.Lanczos sequence for M=3, N=10

Figure 1 :
Figure 1: Lanczos sequence for different configurations of M, and N. Fig.(1a) shows the Lanczos sequence for N=3, and M=3, Fig.(1b) for N=5, and M=3, and Fig.(1c) for N=7, and M=3 respectively for the values of Λ = 0.1, 2.5, 6.6.Fig.(1d) plots the Lanczos coefficients b n for Λ = 2.5, N=10, and M=3.For all the cases, the Lanczos coefficients b n initially grows linearly with n up to n ∼ ln D, and then slowly decrease to zero with a slope of order ∼ − 1 D 2 ∼ − 1 K .The Lanczos sequence {b n } smoothens as we increase the system size.

n Nature of Lanczos coefficients b n 1 1 Kn ∼ D 2 b
<< n < ln D b n grows linearly with n: b n ∼ c n ln D << n < D 2 a steady fall-off of b n with n having slope ∼ n falls off to zero

2 Figure 2 :
Figure 2:The above figure plots the variance σ 2 as a function of the interaction parameter Λ for (N, M ) = (3, 3),(3,5).It can be seen from the plot that as Λ becomes finite starting from Λ = 0, the variance decreases and becomes close to zero for Λ = 1.This behaviour of the variance points towards the transition from integrable to chaotic.

Figure 4 :
Figure 4: Plot of K-complexity at all time scales for Λ = 0.1, 2.5 with N = 5, and M = 3.The linear growth and the saturation can be observed at times of O(D), and at O(D 2 ) respectively.For N = 5, and M = 3, D = 441 and K = 421.The saturation value of C k (t) is K 2 = 220 is shown by the dashed line in the second figure.

bnFigure 5 :
Figure 5: Lanczos sequence Fig.(5a) and K-complexity Fig.(5b) for a different choice of initial operator (A.1).Fig.(5a) shows the Lanczos sequence for N = 5, and M = 3.One can easily see and relate to the previous results and conclude that the Lanczos coefficients b n initially grow linearly with n up to n ∼ ln D, and then slowly decrease to zero with a slope of order ∼ − 1 D 2 ∼ − 1 K .Fig.(5b) plots K-complexity at all time scales for Λ = 2.5 with N = 5, and M = 3.The linear growth and the saturation can be observed at times of O(D), and at O(D 2 ) respectively.For N = 5, and M = 3, D = 441 and K ∼ 421.

Figure 6 : 3 .
Figure 6: Plot of K-complexity at all different scales for Λ = 0.1 with N = 5, and M = 3. Numerical fits are included at both the early time plot in Fig.(6a) near the scrambling time t s = 0.48 and at times ∼ O(D) in Fig.(6b) respectively.Fig.(6a) shows the exponential growth of K-complexity at early times (t ∼ log(lnD)).The non-linear fit function is C K (t) = 0.0379053 exp (0.49635 t) with r 2 = 0.989944 denoted by red dashed lines.Fig.(6b) shows the linear growth of K-complexity at times of O(D).The linear fit function is given by C K (t) = −2.92944+ 1.85039 t with r 2 = 0.992016 denoted by red line.

Table 1 :
The table summarises different behaviors of Lanczos coefficients b n with n for a chaotic system having the dimension of Krylov basis K.

t Nature of Krylov complexity C K
(t) 0 ≤ t ≤ log (ln D) C K (t) grows exponentially in time t : C K (t) ∼ exp c t

Table 2 :
The table distinguishes different behaviors of Krylov complexity C K (t) with t for a chaotic system.