Real-time dynamics of 1D and 2D bosonic quantum matter deep in the many-body localized phase

Recent experiments in quantum simulators have provided evidence for the Many-Body Localized (MBL) phase in 1D and 2D bosonic quantum matter. The theoretical study of such bosonic MBL, however, is a daunting task due to the unbounded nature of its Hilbert space. In this work, we introduce a method to compute the long-time real-time evolution of 1D and 2D bosonic systems in an MBL phase at strong disorder and weak interactions. We focus on local dynamical indicators that are able to distinguish an MBL phase from an Anderson localized one. In particular, we consider the temporal fluctuations of local observables, the spatiotemporal behavior of two-time correlators and Out-Of-Time-Correlators (OTOCs). We show that these few-body observables can be computed with a computational effort that depends only polynomially on system size but is independent of the target time, by extending a recently proposed numerical method [Phys. Rev. B 99, 241114 (2019)] to mixed states and bosons. Our method also allows us to surrogate our numerical study with analytical considerations of the time-dependent behavior of the studied quantities.


I. INTRODUCTION
Many-body localization (MBL) generalizes the concept of Anderson Localization (AL) [1] to the interacting regime and has emerged as a novel paradigm for ergodicity breaking of generic many-body systems subject to strong disorder [2][3][4][5].
In this work, we study the out-of-equilibrium quantum dynamics of bosonic systems deep in an MBL phase in one (1D) and two dimensions (2D). We formulate a method to compute the MBL dynamics of bosons in a controlled approximate fashion by extending a recently proposed method of Ref. 26 to mixed states and bosonic systems. In particular, we focus on local dynamical indicators which are able to distinguish an AL phase from an MBL phase. These dynamical indicators range from temporal fluctuations of local observables to two-point correlation functions, some of which may be used in a cold atom experimental setup, since they involve only few-point correlation functions [20,[27][28][29], and out-oftime-ordered correlators (OTOCs) describing information scrambling in quantum many-body systems [30,31]. We find that bosonic MBL systems exhibit a logarithmic light-cone for information propagation, see Fig. 1 (a,b) for the OTOC in 1D and 2D, respectively. * sunwoo-kim@outlook.com In an MBL phase, particles, despite their interactions, are localized in space, which imposes strong constraints on the dynamics of the system. As a consequence, particle and energy transport are absent and the system retains local information about its past even for asymptotically long times, as opposed to conventional thermalizing systems [15][16][17][18][19][20]. However, for fermionic and spin systems it is by now understood that interactions induce a dephasing mechanism, which allows entanglement and quantum correlations to spread during the dynamics, even though transport is absent [15,32,33]. The entanglement spreading combined with the absence of transport is best understood in terms of an emergent form of integrability and the existence of an extensive number of (quasi-) local integrals of motion (LIOMs) which are coupled by exponentially decaying interactions [34][35][36][37]. As a result, the time evolution produces dephasing of the LIOMs and therefore a logarithmically slow spread of entanglement. On the other hand, due to the local nature of the LIOMs, transport is strongly hindered and charge relaxation is forbidden [26,32,33,35]. In this work we show that this dephasing mechanism also applies to bosonic MBL matter leading to similar phenomenology as has been found for fermions and spins.
The paper is structured as follows. In Sec. II, we describe the method in detail. After detailing the class of models considered, in Sec. II A, we describe the -bit phenomenology of MBL systems, after which we show how to construct an effective -bit Hamiltonian in Sec. II B. Then, in Sec. II C we discuss methods to compute few-body observables, noting that special Gaussian initial states can be chosen in order to efficiently compute quantities in the case of bosons. We also discuss the range of validity of our method in Sec. II D. Then, in Sec. III, we apply the method to study the disordered Bose-Hubbard model focusing on the following local dynamical indicators. We first study the temporal fluctuations in Sec. III A. Then, we discuss how the two-time commutator is an indicator for quantum information spreading, and study its fluctuation, a two-point correlator, in Sec. III B, as well as its second moment, the OTOC, in Sec. III C. Finally in Sec. III D, we study the dependency of these observables on the occupation amplitude, a property unique to bosonic systems.

II. METHODS
In this Section, we describe the model and the method used through the work. A special emphasis will be given on the numerical method used to perform the dynamics and the method's assumptions and limitations.
Throughout the section, we will consider the following general class of lattice models, whereâ † i (â i ) are the creation (annihilation) bosonic/fermionic operators at site i,n i =â † iâ i . and · , · represent pairs of nearest-neighboring lattice sites. The non-interacting HamiltonianĤ 0 is determined by the single-particle Hamiltonian A ij = A ji , which has single particle eigenfunctions {φ l (i)} l with energies {λ l } l . In the presence of sufficiently strong disorder these are Anderson localized, i.e., φ l (i) ∼ e −|i−i l |/ξ loc , where i l is the center of the l th eigenstate. ξ loc (W ) > 0 denotes the localization length and is believed to be dependent on the disorder strength W such that ξ loc 0 as W → ∞. Thus, H 0 can be diagonalized using Anderson creation (annihilation) operators, defined A.
-bit phenomenology In the limit of weak interactions and strong disorder (U W ), it is believed that the system is in the MBL phase and therefore can be described by an extensive number of local integrals of motion [13,38,39], {Î l } l , aŝ In particular, for short-range interactions such as the case in Eq. (1), the interacting LIOMs are expected to be exponentially localized, that is, whereÂ i is a generic observable with finite support realspace around site i, N is the dimension of the Hilbert space, and ξ loc (W, U ) > 0 is the interacting localization length. In the case of bosons, we can account for the unboundedness of the Hilbert space by replac- , where tr N denotes the trace over number of particles up to N . As a consequence, the couplings J (2) l,l ∼ 1 N tr[Î lÎl Ĥ ] ∼ e −|i l −i l |/ξ loc will also be exponentially localized. Importantly, these exponentially decaying interactions imply a non-trivial propagation of information through the system [40].
It is also argued that the interacting eigenstates (U = 0) are adiabatically connected to the non-interacting ones [35,36]. For example, in Ref. 36, it has been analytically shown that at strong disorder and under some mild assumptions on the energy level-statistics, for a random spin chain in 1D, the system is diagonalized via a sequence of local unitary operations. Crucially, the sequence of local unitary operations connect adiabatically the non-interacting eigenstates with the MBL ones, meaning that the exact eigenstates can be approximated at any target precision using a finite sequence of local local unitary operations controlled by the perturbative parameter ∼ U/W .

B. Constructing the effective -bit Hamiltonian
Therefore, as a first approximation, in the limit of strong disorder and weak-interactions, we can approximate the interacting LIOMs {Î l } l with non-interacting ones {Î l } l , as follows. First, we write the Hamiltonian in terms of Anderson operators, i.e.
We approximate the full Hamiltonian by only keeping terms ofĤ int that commute with the non-interacting integrals of motion, retrieving the effective -bit Hamiltonian:Ĥ where l is the renormalized on-site energy, B l,m = B llmm + sB lmlm is the effective interaction term, and s = 1 for bosons and −1 for fermions [26,41,42]. In Eq. (8), we discarded off-diagonal terms of B lmnk . We are justified to do this as the B lmnk term involves the overlap of exponentially localized wavefunctions, and so any off-diagonal elements will be typically much smaller than the on-diagonal ones. We illustrate this below with an argument in 1D, although it is easily generalizable to 2D. Noting ξ loc ∼ 1/ ln W for W 1, we can expand φ l (x) in powers of W −1 : From Eq. (9), we can see that if two indices are different in B lmnk , we gain a factor of W −1 . Therefore the ondiagonal terms are of order ∼ U , whilst the off-diagonal terms are either of order U/W 2 if two indices are same but rest are distinct or of order U/W 4 if all indices are distinct.
The resulting dynamics of Eq. (8) is that of first-order perturbation theory without considering corrections to the wavefunction. Interestingly, this can be seen as a version of Poincaré-Linstedt perturbation theory, in that only resonant corrections to the energy are considered. Unlike conventional perturbation theory, the Poincaré-Linstedt method generates approximate solutions are accurate for all times [43]. By disregarding the changes in the wavefunction, we avoid secular terms in the energies that would diverge over time, see Appendix .

C. Efficiently calculating observables
With the effective -bit HamiltonianĤ -bit , any fewparticle observable can be obtained from expectation values given by time evolution via a collection of effective quadratic Hamiltonians {F lm } l,m which we will define below. We note, however, that the underlying wavefunction is not a Slater determinant, but a much more complicated object. We will first consider quadratic ob-servablesÔ = l,m O l,mη † lη m , whose expectation value can be obtained if all quadratic matrix elements η † lη m t , or equivalently, η mη † l t are known, then generalize to higher order observables. For completeness, we provide the methods for both fermions and bosons.

Product initial states
First, consider an initial state that is a product state in bare-creation operators, It can be shown via elementary commutation relations that for the effective -bit Hamiltonian Eq. (8), where ∆ lm = l − m ,F lm = k F lm kÎ k is the effective quadratic Hamiltonian, whose matrix element is F lm k =B lk −B mk andB lm = B lm + B ml . Let us define U † t = e +itF lm . Using the property U † t U t = 1 and inserting identities between the creation/annhilation operators, we find that where for any operatorÔ,Ô t = e +itF lmÔ e −itF lm . Because the time evolution is via a quadratic Hamiltonian, and the expectation is on the zero-particle state |0 , which is Gaussian for both fermions and bosons, Eq. (11) can be evaluated using the following version of Wick's theorem, which we state for completeness. Wick's theorem. Let the bracket · be an expectation value on a Gaussian state, and ..,N be sets of annhilation and creation operators respectively, where each operator can be in different basis, and also be time-evolved via a quadratic Hamiltonian. Then expectation value of the product of these operators on the Gaussian state can be evaluated to be This method of calculating observables for product states generalizes easily for higher-order expectation values, since we can always put them in a form where all creation operators are on the right, and the annhilation operators are on the left, after which we can apply Wick's theorem.

Mixed initial states
For fermions, expectation values for product initial states can be calculated via (Slater) determinants, which can be done in polynomial time. For bosons, however, this can only be be done using permanents, which scale exponentially with system size. We now show that if we instead consider mixed 'thermal' Gaussian initial states, which are also more relevant in an experimental setting, where pure states cannot be prepared, we can compute the expectation values in polynomial time.
We consider the following initial density matrix, whereQ = i q ini some quadratic Hamiltonian and Z(Q) = tr(e −Q ) is the partition function. In this case, expectation values can be calculated by a trace of the product ofˆ 0 and the desired operator. First note that again using elementary commutation relations, η † l η m t = e it(∆ lm +B ll −Bmm−B lm ) e itF lm tη † lη m . Thus, the expectation value for initial state Eq. (13) can be found to be [44,45] where s = +1 for bosons, and −1 for fermions, the dash operator is defined to rotate a matrix into the -bit basis, so that for any matrix M , Eq. (14) can be derived as follows. First, note that Because bothQ andF lm are quadratic operators, via the Baker-Campbell-Hausdorff formula [46], there exists some quadratic operatorĈ = k C k#k , where# k = γ † kγ k , andη l = a [V ] laγa for some matrix V . Therefore (16) Since the trace vanishes whenever a = b, the trace is a standard bosonic or fermionic counting problem, evaluated to be where C = diag({C k } k ). Substituting and using the summation over V matrices to transform back into the original basis, we find that where Q = diag({q i } i ). Similarly, The partition function can be evaluated to be Z(Q) = det 1 − se −Q −s . Finally, we can substitute the definition of f to find Eq. (14).
Since only a determinant is required, we note that for the initial state in Eq. (13), the computational power required scales polynomially with L and is independent of the targeted time. Similarly, for higher-order expectation values, we can simply represent the expectation in a form such that all the creation operators to the left, and the annhilation operators to the right. Since the expectation is over a Gaussian state time-evolved by an effective quadratic Hamiltonian, we can use Wick's theorem to simplify into quadratic expectation values, then use Eq. 14 to calculate the result in polynomial time.
By substituting t = 0 into Eq. (14), we can see that we can control the initial expectation value of occupation, f i , by appropriately choosing the quadratic Hamiltonian Q as The initial density matrix is thenˆ 0 = {n} w(n) |n n|, where {n} is the set of all possible real-space occupations, and the weights of each occupation is w(n) = i w(n i ), which can be separated into weight of each site as w(n i ) = f n i i /(1 + sf i ) n i +s , which in the case of bosons, is exponentially decaying with occupation, with the rate controlled by f i .

D. Range of validity in the case of bosons
Having introduced our method, let us now take the chance to explore its range of validity.
With this method we can compute observables for bosons in regimes that are far beyond the capabilities of other computational methods such as exact diagonalization or tensor network methods in terms of time t, system size L or number of particles N . Consequently, we cannot directly compare with such an exact reference. Still, our method is expected to be accurate and controlled as we now aim to argue.
In the seminal work of Basko, Aleiner and Altshuler [5], and confirmed by others [36,37] later, it has been shown that perturbation theory starting from the Anderson eigenstates is stable as long as interactions are sufficiently weak. This stability is deeply rooted in the definition of MBL, which implies that integrals of motions are adiabatically connected to the non-interacting ones. In particular, this stability is a consequence of the existence of the spatioenergetic anti-correlations between eigenstates in the single-particle Anderson model. Singleparticle eigenstates that are close together spatially are in general very different in terms of energy. As a consequence, even if denominators of energy differences are small in perturbation theory, the full perturbative corrections are generally nevertheless suppressed exponentially due to the exponentially small overlap of the tails of the involved single-particle wave functions.
This adiabatic continuity to the noninteracting limit forms the basis for our presented method. On the one hand this allows us to choose as the integrals of motion the noninteracting ones to leading order. On the other hand, the off-diagonal interaction terms can, again to leading order, be discarded as they introduce only weak perturbative corrections, see also the case of fermions [26]. Overall, this makes our approach controlled in the limit of strong disorder provided interactions are sufficiently weak. Let us just note, however, that bosons are always expected to become ergodic at sufficiently high energies [39], which naturally limits the occupation numbers, i.e., the density of bosons, in our considered systems.
Concerning the dynamics addressed in this work it is, however, also a natural question up to which time scale the solution can be considered accurate. For fermions it has been shown already that the approach yields controlled results for all times by comparing to an exact reference 26. For bosons we don't have such an exact reference at hand. However, it is nevertheless possible to estimate the temporal validity. In the end our approach becomes uncontrolled whenever resonant processes, beyond the ones already taken care of, contribute significantly. From the simple perturbative analysis in Sec. II.B we have seen already that such resonances don't appear in low-order perturbation theory. Of course, this a priori doesn't exclude significant resonances at higher orders, although it has been argued in a previous work that an MBL phase is supposed to exist at sufficiently low energy densities [39]. Still, even if resonances occur at higher orders in perturbation theory these would only contribute at longer times so that up to that point our description would still be accurate and controlled.
Let us finally emphasize again that in this work we are aiming at the localization dynamics of a 2D bosonic system at strong disorder. It is important to point out that the existence of a genuine MBL phase in 2D is not well established as compared to the 1D case. On the one hand, recent experiments and numerical simulations have shown the stability of an MBL phase in 2D systems [6,[47][48][49][50][51]. On the other hand, general considerations which take into account rare non-perturbative effects have laid out the possibility that an MBL phase in dimension higher than one might not exist asymptotically at large length and time scales [52][53][54][55]. Importantly, these non-perturbative effects evoke the existence of entropic ergodic bubbles, which destabilise the MBL phase at extremely long times [52][53][54][55][56][57]. As a result, even if MBL in 2D might finally turn out to be unstable, we expect that our approach based on the existence of the LIOMs at strong disorder and weak interactions will still be accurate and controlled for fairly long time scales before the instability might kick in.

III. BOSONIC MANY-BODY LOCALIZATION
Having elucidated the main theory of our method, we now apply it to a specific case, namely to a disordered Bose-Hubbard model in 1D and 2D, We consider an initial mixed state, which is close to a charge-density wave state in the sense of f (x=2k,y) = 1, f (x=2k+1,y) = 0 in Eq. (14). We study the system both in a 1D lattice of size L and a rectangular lattice (2D) of size S = L × L/2. Both cases are subject to open boundary conditions. Moreover, we focus our study in the strong disorder and weak interactions (U = 0.1) regime. Here, we expect the aforementioned method to provide an accurate description over several orders of magnitude in time [26].

A. Temporal fluctuations
One of the crucial difference between an MBL phase and an Anderson localized one is the fact the system relaxes in the long-time limit, even though it does not thermalize. This equilibration is manifest in the expectation values of local observables, which relax to stationary values at long times in the thermodynamic limit, and show decaying temporal fluctuations. To explore this difference, in this subsection, we study the temporal fluctuations of the disordered Bose-Hubbard Hamiltonian.
For concreteness, We quantify the temporal fluctuation in the system with the following observable, where · tav indicates the long-time average, x denotes the spatial site index in the x-direction in both 1D and 2D, and y = 0 for 2D. It has been argued [26,27], that the temporal fluctuations decay algebraically in time in an MBL phase in 1D, ∆n 2 (t) ∼ t −α . In particular, it has been predicted that the exponent α is proportional to the localization length (α ∝ ξ loc ).
Here we apply our method to numerically study the decay of time-fluctuations of local observable in 1D and 2D. We confirm the algebraic decaying in 1D and we refine this result in 2D by showing that it may relax as ∆n 2 (t) ∼ e −β ln 2 (t) .
In particular, we present an analytical argument using the effective model and propose a general asymptotic form for the decay of the temporal fluctuation in any dimensions d, ∆n 2 (t) ∼ e −β d ln d (t) . This result, as well as the further analytic predictions that we will present, shows that our method allows us to not only numerically compute observables at regimes inaccesible to other methods, but also allows us to analytically explain the behaviour of these observables. Figure 2 (a), (b) shows ∆n 2 (t) for the 1D and the 2D case, respectively, where · denotes disorder averaging. We can see that in both 1D and 2D, in the interacting case, ∆n 2 (t) decays with time, unlike the non-interacting case U = 0 (dashed lines). The decay for the 1D case is consistent with an algebraic decay ∆n 2 (t) ∼ t −α , that is the dynamical range of time for which a polynomial decay is visible increases with system size (see Fig. 2 (a)). In the inset of Fig. 2 (a), we can also see that α ∝ 1/ ln W , providing evidence that α ∝ ξ loc . In the 2D case, we observe deviations from a pure algebraic decay, and see that the relaxation is consistent with ∆n 2 (t) ∼ e −β ln 2 (t) .
We now present an analytical argument to support the found decay for ∆n 2 (t). First, we assume that the qualitative quench behaviour will not change given similar initial conditions. So, we choose the initial state that corresponds to taking the trace over the (0, 1) site occupation sector, On this state, we find that where˜ q = q + B qq + k =l,mB qk /2. Since F lm k ∼ U e −min(|i l −i k |,|im−i k |)/ξ loc , factors in the cosine product such that i k is further away from i l or i m than ξ loc will be nearly 1, we can restrict the limits product to over k that meets the condition 1 U e −min(|i l −i k |,|im−i k |)/ξ loc t/2, and approximate the factors as some constant less than one. The number of such factors is approximately the volume of two d-balls of radius ξ loc ln(U t). Therefore, there are ∼ ξ d loc ln d (U t) such factors, which points towards the quasi-polynomial decay, as claimed [22,26].

B. Two-point correlators
In this subsection, we study the behavior of two-point correlator function, the fluctuation of the two-time commutator. The two-time commutator of local observables, such as [n x (t),n 0 ], is a key object that describes spatiotemporal propagation of quantum correlations, and can be used to study the spatial growth of local operator such asn x (t), initially localized at site x. This object has also been studied extensively in the context of Lieb-Robinson bounds [58]. To get a feeling for the typical behaviour of its expectation value, we again evaluate on the initial state Eq. (23), and find that [n x (t),n 0 ] ∼ l,m,n,k M lmnk , where the typical matrix element looks like where˜ q = q + B qq + k =l,m,n,kB qk /2. The difference in complex exponential factor suggests that after disorder averaging, the quantity will average to zero. Hence, we can study its fluctuation, a two-point correlator G(x, t): From cosine factor of the typical matrix element, we can use a similar argument as that of ∆n 2 (t) to argue that G(x, t) too will show a quasi-polynomial decay at moderate times.
In the 1D case, Fig. 1 (c) shows a color plot of G(x, t), from which we can deduce its qualitative properties. In the 1D case, Fig. 1 shows that in the case of the AL phase, we see an initial spread of fluctuations up to localisation length, followed by a static phase, which persists over many decades. Meanwhile, the MBL phase has a different qualitative behaviour; at small times, G(x, t) spreads and grows, but unexpectedly, starts to decay at moderate times [26]. The quantitative growth of G(x, t) is shown in Fig. 3 (a). In the AL case, we see that there is an initial growth phase followed by saturation. For the MBL case, however, there exists a time t * , around which the initial growth phase turns into what appears to be a polynomial decay, as predicted by the analytic argument. We expect that t * ∼ 1/U , on dimensional grounds and the fact that the peak is due to interactions, since it does not exist in the Anderson-localized case. This argument is supported by the inset of Fig. 3 (a). Additionally, the data collapse from adding a factor of |x| from G(x, t) indicates that the magnitude of G(x, t) decays exponentially in space, i.e. G(x, t) ∼ e −|x|/ e −C d ξ d loc ln d (U t) . We expect that ∼ ξ loc , as ξ loc is the only relevant length scale of the system. This is again supported by the inset of Fig. 3 (b). Finally in Fig. 1 (d) and Fig. 3 (b) show the same phenomena, but for the 2D case, except for the case of Fig. 3 (b), where the decay is consistent with a quasi-polynomial decay, i.e. ∼ ln 2 (t) for ln(t) 0.

C. Out-of-time-order correlators
In this subsection, we inspect the behavior of the OTOC in 1D and 2D. The OTOC has been widely studied in the context of MBL [59,60] and is thought to describe quantum information scrambling in many-body systems, and is known to diagnose the onset of chaos the semi-classical limit [31,61]. It is also the second moment of the two-time commutator discussed in the previous subsection and is defined as The OTOC, which is 8th-order in creation and annhilation operators, has been numerically simulated for spin chains and fermions, but never for a bosonic system at the time of writing to the best of our knowledge. Our simulation of the OTOC shows that our method can be used to simulate high-order correlation functions. Fig.  1 (a), (b) shows a color plot of C(x, t) for 1D and 2D respectively. We see that for the AL case, the OTOC spreads only up to the localisation length then becomes stationary -correlations no longer spread over time. On the other hand, the MBL phase shows stark difference to the AL phase, where we can clearly see a logarithmic light-cone that persists even at very large times, indicating that quantum correlations continue to spread at constant 'log-velocity' c with ln(t). In Fig. 4 (a), we see that the OTOC first has a growth phase at a power law of ∼ t 2 , then then saturates to a constant value at long times. This agrees with our analytic prediction, which we will present below, as well as those of spin chains [30,59], providing evidence that OTOCs in bosonic systems share similar features to those of spin chains. The data collapse by shifting of time via log-velocity c, as shown in Fig. 4 (c), confirms that the spread is indeed linear in ln(t). In the inset of Fig. 4 (c), we study how the logvelocity c scales as a function of disorder strength W . We see that the log-velocity scales as 1/ ln W , which, since ξ loc ∼ 1/ ln W for strong disorder, means that increasing localization length increases the log-velocity of the light cone. Much of the same observations can be seen in 2D from Figs. 1 (b), 4 (b) and 4 (d). We now provide arguments for the observed spatiotemporal behaviour of the OTOCs. In a non-interacting localized phase, we can expect that after an initial spread, C(x, t) will have support for |x| ξ loc . This is because the non-interacting operatorsη † l = i φ l (i)ĉ † i , which the system evolves with, only have support on sites ξ loc around i l , which is what we saw. On the other hand, in the interacting phase, we expect thatn x (t) will have significant support at site 0 only when t ∼ 1/B x0 . Since B lm ∼ U e −|i l −im|/ξ loc , it suggests non-zero C(x, t) for |x| ξ loc ln(U t). Thus, we expect a logarithmic light-cone ∼ ln(t), with the log-velocity scaling with localization length c ∼ ξ loc , as observed. We suspect that the reason this argument does not hold for G(x, t) is because the decay process is the faster, more dominant process, compared to the slow spread of correlation. We can also argue the early-time growth of the OTOC should go as ∼ t 2 , as follows. First, we assume that the qualitative behaviour should not depend on the choice of the local operators, as long as they are not correlated at t = 0. So, we choose the commutator [Î l (t),Î 0 ], whereÎ l = m,k∈{l,l+1}η † mηk , such that |i l | ≥ ξ loc . Next, we argue that the early time scaling behaviour can be deduced from one of the terms, since the OTOC is just a sum of them. Choosing one of the terms, it can be shown that where˜ q = q + B qq + k =0,l,l+1B kq /2 +B q0 . Upon expanding this expression for small t, we retrieve the ∼ t 2 growth at early times, consistent with the observed result.

D. Dependency on occupation amplitudes
Now that we have studied properties shared between bosonic, fermionic and spin-chain systems, we now consider the effect of the occupation amplitude f amp on the studied observables {O}, a property unique to bosons. We consider Néel-like initial states with various amplitudes f (x=2k,y) = f amp , f (x=2k+1,y) = 0. First, we study (1 + δA) n = 1 + nδA + O(δ 2 ), to find that quadratic observables should scale as ∼ f amp . For the case f amp 1, we can utilize the semi-classical limit to approximate expectations of occupations as simply f amp . Since ∆n 2 (t) is a square of a quadratic observable, we expect it to scale as ∼ f amp 2 for both limits. For the other two observables G(x, t) and C(x, t), however, we must exercise care as they contain commutators. We can make predictions about these observables in the semi-classical limit f amp 1. Generally, the commutator 'eats up' operators, leaving a commutator of two operators of as the same order as the commutands themselves So, we can approximate the commutator as a random variable of order f amp . Since the mean vanishes, we expect it to be distributed ∼ ±f amp , which means that the absolute value has the expectation value f amp for each disorder realization. Since G(x, t) is the square of the absolute value of the commutator, we predict a linear scaling G max ∼ f amp . On the other hand, the OTOC measures the square of the commutator, and therefore has the expectation value f amp 2 for each disorder realization. Therefore, we can predict a C max ∼ f amp 2 for f amp 1. In Fig. 5 (a), we show the scalings of O max for the three studied observables. We see the scaling for ∆n 2 max is very close to ∼ f amp 2 , as predicted by the analytic arguments. G max is also consistent with a linear scaling, ∼ f amp . We note that surprisingly, ths scaling holds for f amp < 1 as well. Lastly, we see that C max appears to be consistent with interpolation between linear and quadratic behaviour in f amp .
Next, we study how the power-law decay constant α changes with occupation amplitude f amp . To come up with an analytic argument, we choose a similar state as Eq. (23), but instead of a trace up to one particle, we take a trace up to f amp particles. For this state, instead of the cosine factor that we see in Eq. (24), which arises from summing e −itF lm k n for n = 0, 1, we instead get a truncated geometric series for n = 0, . . . , f amp . The magnitude of the equivalent factor is ∼ |sin((f amp +1)F lm k t/2)/ sin(F lm k t/2)|, whose maximum value scales as ∼ f amp . Therefore, we can approximate factors inside of the d-ball as a constant proportional to f amp , Cf amp . This results in the time-dependent power law to scale as α ∼ ln(Cf amp ) ξ loc . However, in the limit f amp → 0, we expect α → 0. This is because dynamics are caused by interactions, which is suppressed if there are no particles. This argument is numerically supported in Fig. 5 (b), where we see that plot of exp(α) vs. f amp results in a straight line, for both ∆n 2 (t) and G(x, t).

IV. CONCLUSION
In this work, we have formulated a method to solve for the dynamics in 1D and 2D bosonic quantum matter in a deep MBL phase at strong disorder and weak interactions. This method represents an extension of the approach described in Ref. [26] and yields an algorithm requiring resources depending polynomially on system sizes but being independent of the targeted time. We find that this allows us to access regimes far beyond the capabilities of other methods such as exact-diagonalization or tensor networks. Using this method, we have provided evidence that many of the MBL indicators, such as the OTOC and temporal fluctuations, which were studied and well established for fermionic and spin-chain systems, also apply to bosonic MBL matter both numerically and with analytic arguments. Our approach can also be directly applied to other nonergodic settings such as Stark MBL [62,63]. We have checked that the resulting physical properties very much resemble the disordered bosonic system studied in this work.
While we have been concentrating mainly on few-body observables up to OTOCs, it remains an open question to solve for the dynamics of more complicated quantities involving also higher-order correlation functions such as entanglement measures. Computing these with our method appears challenging as entanglement measures for interacting theories cannot be reduced to single-particle observables in general, especially when targeting entanglement of larger subsystems. A potential route towards accessing such quantities could be to use an approach based on correlations between random measurements [64,65], which has also been used recently in an experiment to measure second-order Renyi entropies [66].
The present work might be considered as a first step towards exploring the dynamics of bosonic MBL matter. In the future it might be interesting to increase the accuracy of the approach. This might be achieved by accounting for the dressing of the LIOMs, which currently are taken to be the Anderson orbitals, due to interactions. A crucial next step might also be to take into account the offdiagonal interaction terms. In principle, they might be integrated out using a Schrieffer-Wolff transformation at the expense of generating higher-order diagonal terms. The resulting dynamics, however, cannot be captured anymore by the presented free-boson techniques as now the few-body correlation functions cannot be evaluated by means of Wick's theorem.
In a more general context, our approach shares similarities with the recent ideas to use classical [67,68] or artificial neural network wave functions [69][70][71] for the solution of Schrödinger's equation. Here, our method can be interpreted as a perturbative construction of such networks upon explicitly designing the network structure and computing the respective weights directly without any training or optimization [26]. Pushing our approach further to include for instance the off-diagonal interaction terms might be done by taking this connection to classical or artificial neural network wave functions and to use a time-dependent variational principle to optimize our network parameters by means of stochastic Monte-Carlo techniques. In the end it might be possible to also get closer to the experimental regimes of larger interactions and potentially towards explaining or developing a theoretical description of the pioneering experiment on bosonic MBL matter [6]. ACKNOWLEDGEMENT We would like to thank M. Knap for enlightening discussions. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 853443), and M. H. further acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) via the Gottfried Wilhelm Leibniz Prize program.
Appendix: Poincaré-Lindstedt perturbation theory In the non-interacting limit U = 0, the Hamiltonian of Eq. (6) is nothing but a collection of decoupled harmonic oscillators. The presence of interactions introduces a coupling between these oscillators in a non-linear fashion. In our approximation in Eq. (8), we only consider the diagonal contribution of the interaction. This takes into account the frequency shift of the oscillators and therefore resembles Poincaré-Lindstedt perturbation theory, which is a powerful perturbative approach for solving harmonic oscillators in classical systems.
Let us therefore take the chance to illustrate the working principle of our approximation using this analogy. For concreteness, let's consider the so-called Duffing equation, the simplest non-linear perturbation to the simple harmonic oscillator: where 0 < 1 accounts for a weak nonlinearity. The Poincaré-Lindsted method is a version of perturbation theory, where, by considering the shift in the frequency of the unperturbed oscillatory solution, secular terms are systematically removed, resulting in a oscillatory solution that is accurate for all times.
(A.3) We see that from the first order in , a secular term ∼ t sin(t) arises, which is not bounded in time.
In Poincaré-Lindstedt perturbation theory, on the other hand, in addition to the amplitude, we also consider the shift in frequency, as τ = ωt, where ω = 1+ ω 1 +· · · . By appropriate choice of ω 1 , the term proportional to t sin(t) can be removed, resulting in the solution x PL (t) = cos 1 + 3 8 t (A. 4) for the case of the Duffing equation. Within our approach we go one step further. We only consider frequency shifts and neglect any other perturbative contribution. For the above Duffing oscillator this amounts to neglecting x 1 (t) (and higher-order corrections), resulting in the approximation x Ours (t) = cos 1 + 3 8 t . (A.5) In Fig. 6, we compare x Ours (t), a simplification of the Poincaré-Lindstedt perturbation theory (blue line), with the conventional perturbation theory (orange line) and the exact solution (magenta dashed line). We see that even though Our solution neglects the term proportional to in the Poincaré-Lindstedt solution, it is still accurate for moderate perturbative strength and stays bounded over time, whereas the conventional perturbation theory breaks down for times t 1/ due to the presence of secular terms.