Localization dynamics in a centrally coupled system
Nathan Ng, Sebastian Wenderoth, Rajagopala Reddy Seelam, Eran Rabani, Hans-Dieter Meyer, Michael Thoss, Michael Kolodrubetz
LLocalization dynamics in a centrally coupled system
Nathan Ng, ∗ Sebastian Wenderoth, ∗ Rajagopala Reddy Seelam, EranRabani,
3, 4, 5
Hans-Dieter Meyer, Michael Thoss,
2, 7 and Michael Kolodrubetz Department of Physics, University of California, Berkeley, CA 94720, USA Institute of Physics, University of Freiburg, Hermann-Herder-Strasse 3, 79104 Freiburg, Germany Department of Chemistry, University of California, Berkeley, CA 94720, USA Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA The Sackler Center for Computational Molecular and Materials Science, Tel Aviv University, Tel Aviv 69978, Israel Theoretische Chemie, Physikalisch-Chemisches Institut,Universit¨at Heidelberg, INF 229, D-69120 Heidelberg, Germany EUCOR Centre for Quantum Science and Quantum Computing,University of Freiburg, Hermann-Herder-Strasse 3, 79104 Freiburg, Germany Department of Physics, The University of Texas at Dallas, Richardson, Texas 75080, USA (Dated: February 2, 2021)In systems where interactions couple a central degree of freedom and a bath, one would expectsignatures of the bath’s phase to be reflected in the dynamics of the central degree of freedom.This has been recently explored in connection with many-body localized baths coupled with acentral qubit or a single cavity mode — systems with growing experimental relevance in variousplatforms. Such models also have an interesting connection with Floquet many-body localizationvia quantizing the external drive, although this has been relatively unexplored. Here we adapt themultilayer multiconfigurational time-dependent Hartree (ML-MCTDH) method, a well-known treetensor network algorithm, to numerically simulate the dynamics of a central degree of freedom,represented by a d -level system (qudit), coupled to a disordered interacting 1D spin bath. ML-MCTDH allows us to reach ≈ lattice sites, a far larger system size than what is feasible withexact diagonalization or kernel polynomial methods. From the intermediate time dynamics, we finda well-defined thermodynamic limit for the qudit dynamics upon appropriate rescaling of the system-bath coupling. The spin system shows similar scaling collapse in the Edward-Anderson spin glassorder parameter or entanglement entropy at relatively short times. At longer time scales, we see slowgrowth of the entanglement, which may arise from dephasing mechanisms in the localized systemor long-range interactions mediated by the central degree of freedom. Similar signs of localizationare shown to appear as well with unscaled system-bath coupling. I. INTRODUCTION
The advent of controllable quantum simulation plat-forms allows for novel explorations of quantum coherentphenomena. Certain such architectures have the advan-tage of using extra degrees of freedom as a way to easilyread out properties of a system [1]. Examples of suchsetups include cavity QED with ultracold atoms [2] andsuperconducting qubit circuits, the latter of which was re-cently used to simulate the many-body localized (MBL)phase in a 10 qubit chain with long-range interactionsmediated by a central resonator [3]. Given that suchplatforms are in their early stages, it is important to ex-plore the interplay of disorder-induced localization andmediated long-range interactions, and how they affectthe dynamics of localization in these systems.If localization exists in these systems, it will naturallybe many-body localization since the spins hybridize withthe central degree of freedom to give non-trivial interac-tions. Rigorous results on MBL have already been es-tablished in one dimensional systems with short rangedinteractions [4]. In such a setting, it is a stable phase ∗ These authors contributed equally to this work. of matter, with respect to adding short range perturba-tions, that can coexist with other types of order [5, 6].While strong disorder enables localization, it cannot pre-vent thermalization if interactions are long-ranged, de-caying slower than r − D , where D is the spatial dimen-sion [7, 8]. Even the MBL phase with short-ranged in-teractions is fragile. It is destroyed upon coupling toa continuum of bath modes [9] which, intuitively, canprovide arbitrary amounts of energy and allow the sys-tem to transition between eigenstates of vastly differentcharacter. One sees then that there are two ingredientsto this delocalization mechanism: a continuum of ener-gies of large enough bandwidth, and hybridization dueto effective infinite-ranged interactions mediated by thenon-Markovian bath.In fact, for a specific type of memoryless bath, noner-godicity does survive. This is the case of Floquet MBL,in which an MBL system is subjected to an external peri-odic drive with frequency Ω modeled as a time-dependentHamiltonian acting on the system [10, 11]. The failure ofthermalization is due to the the inability of the systemto absorb energy in quanta of (cid:126) Ω, which itself is a conse-quence of the discreteness of the energy spectrum. Theexternal drive, however, is not inherently dynamical andthus does not capture the backaction present in a fullyquantum mechanical system. a r X i v : . [ c ond - m a t . d i s - nn ] J a n In this work, we consider the time evolution of such asystem obtained by treating the Floquet drive as a quan-tum degree of freedom. Specifically, we consider a local-ized system globally coupled to a d -level system (qudit)with finite energy spacing, similar to [12]. When thequdit is a two-level system, it was shown that localiza-tion does not survive at any finite coupling [13, 14]. Butwhen it is instead a d > II. MODEL
We consider a simplified model of many-body localiza-tion by coupling a one-dimensional chain of qubits (spins-1 /
2) via global interactions with a central qudit: H = H + Ωˆ τ z + γH (cid:0) ˆ τ + + h.c. (cid:1) , (1) H = L (cid:88) i =1 hξ i σ z i + gσ z i σ z i +1 , H = L (cid:88) i =1 σ x i , where ˆ τ z = (cid:80) dn =1 n | n (cid:105)(cid:104) n | , ˆ τ + = (cid:80) d − n =1 | n + 1 (cid:105)(cid:104) n | , andthe operators H and H act only on the spin subspace.The states | n (cid:105) label the states of the central qudit. Here, h = 1 . g = 1 .
07, Ω = π/ .
8, and ξ i is a random vari-able drawn uniformly from ( − , d → ∞ ), it shows a localization-delocalization transition at a critical coupling γ c (cid:46) . γ either deep in thelocalized phase ( γ < .
2) or deep in the ergodic phase( γ ≈ d = 7, which is large enoughto display Floquet-like behavior but small enough thatthe finite qudit size plays an important role.The spin part of the Hamiltonian, H , is a trivial an-tiferromagnetic Ising chain with longitudinal on-site dis-order. The diagonal nature of H in the z -basis yieldstrivial localization in the eigenstates. This manifests ineigenstates | ψ n (cid:105) as vanishing site-averaged magnetization L − (cid:80) i (cid:104) ψ n | σ z i | ψ n (cid:105) and maximal value of the spin-glassparameter, q = L − (cid:80) i (cid:104) ψ n | σ z i | ψ n (cid:105) = 1 at high energydensities. Values of q ≈ q → τ z onthe qudit, thus selecting a preferred direction for the cen-tral spin. This greatly impacts the ease with which thequdit fluctuates, which in turn can regulate transitionsin the spin states leading to delocalization.Recent studies have examined how localization canpersist in the presence of long-ranged interactions [8, 22]or with central coupling to a single degree of freedomyielding an effective Hamiltonian with long-ranged inter-actions [13, 14]. With the exception of a numerical study[22], these past works have noted that preservation oflocalization in the thermodynamic limit requires increas-ing the disorder strength with increasing system size ordecreasing the strength of central coupling as γ → γ/L .Reducing the coupling strength in this way renders thelong-ranged part of the effective Hamiltonian for the spinchain subextensive. This is also reflected in the dynamicsof the qudit as its transition rate vanishes.On the other hand, the existence of Floquet MBL af-fords a different pathway to the coexistence of localiza-tion and central coupling. In that context, the persis-tence of MBL is not due to a vanishing coupling to theexternal drive, but to a suppression of mixing betweendifferent localized eigenstates of the undriven system.This picture suggests that an effective Hamiltonian foronly the spin degrees of freedom should show localizedbehavior. This is indeed the case, as previous work basedon the high frequency expansion has shown [12]. In thislimit of Ω → ∞ the spins are governed by an effectiveHamiltonian diagonal in the qudit basis, reproducing theeigenenergies modulo an integer multiple of Ω: H eff = H + ( H ) | d (cid:105)(cid:104) d | − | (cid:105)(cid:104) | Ω + O (Ω − ) . (2)At lowest order in Ω − , we see that possible delocaliza-tion is reserved only for states with | (cid:105) or | d (cid:105) , as ( H ) induces all-to-all coupling. Increasing L without increas-ing d , as we do in this paper, means that eigenstatesoccupying | (cid:105) will eventually encroach upon the middleof the spectrum and contribute to the quench dynam-ics we study. This can be seen from the density of stateswhen H is dominant as it follows ρ ( E ) ∝ exp (cid:16) − E ( J √ L ) (cid:17) for energy scale J ∼ O (1), meaning ρ ( E ) will grow widerwith increasing L . An energetically dominant ( H ) termwill both delocalize the eigenstates and deform the Gaus-sian density of states in the thermodynamic limit.We thus assume γ to be small enough such that neitheroutcome occurs, and ask when this picture will naivelybreak down. In such a limit, we can treat the H fieldterm in a mean field fashion for each eigenstate: H eff ≈ H + γ L Ω + (cid:88) i γ Ω (cid:42)(cid:88) j (cid:54) = i σ x j (cid:43) σ x i , where the effective field (cid:68)(cid:80) j (cid:54) = i σ x j (cid:69) in an eigenstate mustbe determined self-consistently. For a typical eigenstate,this field should have value ∼ f ( γ ) √ L , where f ( γ ) mustvanish when γ = 0. This is the case when (cid:80) j (cid:10) σ x j (cid:11) is thesum of L − γ eigenstates are assumed to be perturbatively connectedto a corresponding γ = 0 eigenstate. For this model, wetake the lowest order approximation f ( γ ) ≈ f γ . Withthis assumption [24], the effective transverse field on site i will begin to compete with the longitudinal fields in H when γ (cid:68)(cid:80) j (cid:54) = i σ x j (cid:69) ∼ O ( g, h i ) ∼ O (1). For the highenergy density eigenstates we are interested in, this ef-fective field will inhibit spin glass ordering and the sys-tem should obey the eigenstate thermalization hypothe-sis. Thus, γ ∝ L − / should serve as a rough separatrix d √ Lγ ThermalizingLocalizing0 1 ∞ γ ∗ c . . .
35 III L = 24 L = 48 L = 96 FIG. 1. Schematic phase diagram for the coupled system (1),along with the parameters for which we present numericalresults from ML-MCTDH. The rough phase boundaries aredetermined from numerics and analytical arguments. Thedata are separated into three solid segments – the strong,intermediate and weak couplings from top to bottom. Theangle of the segments comes from fixing the qudit size to d = 7 and scaling the central coupling γ ∝ L − / . For easeof discussion, we group the three coupling regimes as region I(weak and intermediate) and region II (strong). between thermalizing and athermal behaviors. Further-more, couplings that tend to zero faster than L − / willrealize a trivial limit, where the localization comes en-tirely from H . The region where this is argument ex-pected to be most significant is denoted in Fig. 1 througha color gradient starting around d/ √ L ∼
1. Note that,in general models where H includes operators diagonalin the z-basis, we would have f (0) (cid:54) = 0; in this case thescaling is replaced by γ ∼ L − / .Besides scaling the coupling to zero, the all-to-all in-teractions can be avoided by ensuring that eigenstatesoccupying levels | d (cid:105) or | (cid:105) in the qudit do not partici-pate in the dynamics. For quenches starting from themiddle of the many-body spectrum, this condition canbe ensured by keeping the qudit size d sufficiently largecompared to the typical width of the 1D many-body den-sity of states, √ L . Dynamics in this limit should closelyresemble Floquet physics, since the fluctuations produc-ing effective long-ranged interactions will cancel out af-ter accounting for the processes in which the intermedi-ate qudit state changes by +1 or −
1. Away from thislimit, when d/ √ L (cid:46) O (1), the all-to-all interactions areunavoidable. The threshold value of d/ √ L for delocaliza-tion should decrease as the coupling is decreased. Thesearguments are summarized schematically in Fig. 1.There are two important ways to think of this systemand its dynamics: either as a combined many-body sys-tems with localized and delocalized phases, as was donein the previous paragraph, or as a central qudit inter- | Ψ (cid:105) = N (cid:88) j =1 ... N (cid:88) j P =1 A j ,...,j P ( t ) P (cid:89) κ =1 | ϕ ( κ ) j κ ( t ) (cid:105) , | ϕ ( κ ) j κ ( t ) (cid:105) = N (cid:88) i =1 ... N (cid:88) i Q ( κ ) =1 B κ,j κ i ,...,i Q ( κ ) ( t ) Q ( κ ) (cid:89) q =1 | ν ( κ,q ) i q ( t ) (cid:105) ,..., | Ψ iQ N = d bcd S S S N d S S S N N bc S S S N c S S S N N N = d FIG. 2. Expansion of the wave function | Ψ (cid:105) and the first layer single-particle functions | ϕ ( κ ) j κ ( t ) (cid:105) used in the ML-MCTDHapproach (left) and a schematic representation of the tree structure of the wave function (right) . The black dots representsingle-particle functions (SPFs). The red dot represents the qudit degree of freedom and the blue dots represent the spin degreesof freedom. The binary expansion of the spin wave function is symmetric and, thus, we choose the numbers of SPFs within onelayer to be equal. In the example shown, only three spins are grouped together in the lowest layer for better visualization. Inthe calculation, however, groups of up to 12 spins in the lowest layer are used. acting with an unusual, localized, spin bath. From thislatter viewpoint, it will be useful to consider scaling thesystem bath coupling γ ∼ / √ L , since that will be shownto achieve a well-defined thermodynamic ( L → ∞ ) limit.This scaled coupling will be used in the majority of oursimulations, and is covered in more detail in Section IV.For now, we note that γ ∼ L − / scales to zero fasterthan the L − / that we predict is required for MBL.Therefore, at sufficiently late times, we predict MBL withour scaled coupling.In this model we use qudits for numerical simplicitydue to their finite Hilbert spaces. However, our conclu-sions can be easily applied also to the case where thecentral degree of freedom is a single bosonic mode, suchas in cavity QED or superconducting circuits. In thesesetups we expect similar dynamical behaviors when thecentral coupling is appropriately scaled [12]. III. NUMERICAL METHOD
The non-local interaction induced by the centrally cou-pled qudit makes the simulation based on matrix prod-uct operator techniques like time-evolving block decima-tion [25] inefficient. And while alternative approachessuch as the Floquet-Keldysh DMFT [26] exist, they arevalid only in the well-studied Floquet limit in which theinteresting mediated all-to-all couplings are negligible.Thus, we instead employ the Multilayer Multiconfigu-ration Time-Dependent Hartree (ML-MCTDH) method[15, 16, 18, 19, 27] which has been used to study simi-lar systems in the past, e.g. a two-level system coupledto a bath of noninteracting spins [28]. The ML-MCTDHmethod generalizes the original MCTDH method [17, 29–32] for applications to significantly larger systems. TheML-MCTDH approach represents a rigorous variationalbasis-set method, which uses a multiconfiguration ex-pansion of the wave function, employing time-dependent basis functions and a hierarchical multilayer representa-tion. Within this framework the wave function is recur-sively expanded as a superposition of Hartree productsas depicted in Fig. 2. Here, | ϕ ( κ ) j κ ( t ) (cid:105) , | ν ( κ,q ) i q ( t ) (cid:105) , . . . , are the so-called “single-particle functions” (SPFs) forthe first, second, etc. layer and the coefficients A j ,...,j N , B κ,j κ i ,...,i Q ( κ ) are the expansion coefficients of the first, sec-ond, etc. layer. Despite their name, the SPFs describemultiple degrees of freedom, see Fig. 2. The ML-MCTDHequations of motions for the expansion coefficients andthe single-particle functions are obtained by applying theDirac-Frenkel variational principle [15, 33], thus ensur-ing convergence to the solution of the time-dependentSchr¨odinger equation upon increasing the number ofSPFs. In principle, the recursive multilayer expansion,which corresponds to a hierarchical tensor decomposi-tion in the form of a tensor tree network, can be carriedout to an arbitrary number of layers. In practice, themultilayer hierarchy is terminated at a particular levelby expanding the single-particle functions in the deepestlayer in terms of time-independent basis functions.In the present application of the ML-MCTDH method,we separate the qudit wave function and the spinchain wave function in the uppermost layer as depictedschematically in Fig. 2. The wave function of the spinchain is then further expanded in a binary tree (i.e. P = Q = 2) up to the lowest layer, which comprisesblocks of up to 12 spins. Each of the lowest blocks isexpanded in the time-independent local basis of the un-derlying Hilbert space. Regarding the number of SPFsin the first layer, N , it can be shown that N > d leadsto redundant configurations in the expansion [29], andthus, we set N = d in all calculations. The requirednumber of SPFs in the other layers of the expansion ofthe spin-chain wave function was determined by thoroughconvergence tests and depends on the coupling strength γ . In general, fewer SPFs are needed for smaller couplingstrengths. For L=24, two dynamical layers are employed − ∆ Q / . (a) L=12L=24 L=48L=96 − L ( − q E A ) / . (b) − . . . . S A (c) − Ω t/ (2 π ) ∆ Q / . (d) Noninteracting spins ( g = 0)Perturbation theory − Ω t/ (2 π ) L ( − q E A ) / . (e) − Ω t/ (2 π ) S A (f) − Ω t/ (2 π )0 . . q E A − Ω t/ (2 π )0 . . . q E A FIG. 3. Dynamics from an initial “super-Neel” state in the regimes of weak (top row, γ = 0 . / (cid:112) L/
12) and intermediatecoupling (bottom row, γ = 0 . / (cid:112) L/ O (10 ) – O (10 ) disorder realizations, with the shadedbands indicating deviations of ± (left) Variance of qudit occupations ∆ Q = (cid:10) ( τ z ) (cid:11) − (cid:104) τ z (cid:105) , (center) deviation from perfect spin glass order 1 − q EA , and (right) entanglement entropy S A between a contiguous halfof the spin chain and its complement. The observables ∆ Q and 1 − q EA have been appropriately rescaled [24] to show theircoincidence (except a factor of 2) at early times for weak coupling. The dynamics are observed to converge to a single curve(black dotted line) and appear to be consistent with the dynamics without the nearest neighbor Ising coupling (dashed lines)as L → ∞ . and the required number N of SPFs varies from 30 to120 SPFs. For L=48, a three layer scheme is used wherethe number of SPFs in the lowest layer varies from 10to 30 and in the highest layer from 20 to 60 SPFs. ForL=96, four layers are employed with SPFs which varyfrom 10 to 20 in the lowest and from 35 to 50 in thehighest layer. IV. RESULTS FOR SCALED COUPLING
We examine the system at infinite temperature by fo-cusing on states in the middle of the many-body spec-trum, which have energies close to the midpoint betweenthe maximal and minimal energies of the coupled system,( E max and E min ) respectively. We take γ = 0 for t < | ↓↓↑↑ . . . (cid:105) and thequdit occupying its middle state | ( d + 1) / (cid:105) . The cou-pling is switched on instantaneously at t = 0 to a finitevalue. The super-Neel state is on average a zero energyeigenstate of H and has subextensive energy variance,making it a suitable microcanonical probe. Thus, whenthe system is thermalizing and shows ensemble equiv-alence, we expect similar dynamics compared to onesobtained through averaging over random initial productstates, mimicking an infinite temperature canonical en-semble.As there are different dynamical behaviors in ourmodel, we shall organize our discussion around theschematic phase diagram in Fig. 1, similar to the onefirst introduced in [12]. In this first section, we will con- sider scaling the coupling as γ ∼ / √ L , correspondingto the three solid lines in the phase diagram which, fromtop to bottom, will be referred to as the strong, inter-mediate, and weak coupling regimes. The orientation ofthese cuts comes from the 1 / √ L scaling of γ . This isnatural if we think of the qudit as our main object of in-terest, as it gives a well-defined thermodynamic limit forthe qudit when it is coupled to a non-interacting bath,such as in the spin-boson model [34]. This scaling repro-duces the Kac prescription [35] for the all-to-all term inthe effective Hamiltonian ensuring also the existence of athermodynamic limit for the spins. Specifically, we scale γ using the following formula: γ = γ (cid:114) L L (3)where L = 12 throughout for convenience, such that γ = γ at L = 12. γ sets the overall strength of the coupling.We will consider three regimes, indicated by the solidlines in Fig. 1: weak coupling ( γ = 0 . γ = 0 . γ = 0 . A. Weak and intermediate coupling (region I)
The first cases we consider are weak and intermediatecoupling, which are labeled Region I in Fig. 1. Theseare both at sufficiently small γ that we expect MBL forthe largest accessible system sizes, but for intermediatecoupling ( γ = 0 .
3) the system will be near the phase − ∆ Q / . (a) − L ( − q E A ) / . (b) L=8L=12L=16 L=24L=48L=96 − . . . S A (c) − ( p n + p d − n ) / × − (d) | n − n mid | = 1= 2= 310 − γ t/ Ω ∆ Q / . (e) − γ t/ Ω L ( − q E A ) / . (f) − γ t/ Ω S A (g) log(0 . γ Ω t ) + 1 . − γ t/ Ω . . . ( p n + p d − n ) / (h) | n − n mid | = 1= 2= 3 FIG. 4. Same as in Fig. 3, but with time rescaled by the system size-dependent coupling γ . Grey dots are independentcalculations using the kernel polynomial method, aimed to extend the maximum time from Ω t/ (2 π ) ∼ × to Ω t/ (2 π ) ∼ . × . Between the weak (a,b,c,d) and intermediate (e,f,g,h) coupling regimes, there is a qualitative shift in the longtime behavior of both the qudit- and spin-only observables. This data is suggestive of logarithmic growth in the entanglemententropy becoming the dominant characteristic after t ∼ /γ . (d,h) Occupations p n of the qudit levels, symmetrized aroundthe middle level | n mid (cid:105) = | (cid:105) . While both the weak (d) and intermediate (h) coupling regimes have most of their populationsconcentrated the initial occupied level, | (cid:105) , the latter case has a much greater fraction of the total population in the extremesof the qudit’s states. The values of p n for | n − n mid | = 3 in the left panel are too small ( ∼ O (10 − − − )) for the scale. transition for small L . Three observables – the quditvariance ∆ Q (Eq. 4), the spin glass order parameter q EA (Eq. 5), and the entanglement entropy of the half chain S A (Eq. 6) – are plotted in Figs. 3 and 4, which corre-spond to identical data with different scaling of the timeaxis. The origin of this scaling will be clarified shortly.Note first that, by preparing both the qudit and thespins in highly excited states, one would normally expectthe system to relax quickly to a featureless “infinite tem-perature” equilibrium. That is, all internal levels of thequdit should be equally occupied, and the spins should beparamagnetic and translationally invariant. This is nottrue for the disordered system we study, as the numericsdemonstrate in Fig. 3: for sufficiently small coupling, thesystem shows localization in both the qudit and its thesurrounding spins. The former is signaled by the varianceof the qudit occupations∆ Q ≡ (cid:10) (ˆ τ z ) (cid:11) − (cid:104) ˆ τ z (cid:105) , (4)which saturates to a quantity far below that of the uni-form limit, ∆ Q = ( d − /
12 = 4. Furthermore, thedifferent system sizes exhibit scaling collapse of ∆ Q upto a time scale t ∼ /γ . This is a property of the scaled γ , as it implies that the spin chain acts as a bath for thequdit with a well-defined thermodynamic limit. Morespecifically, it can be shown that the considered modelwith scaled coupling γ ∝ L − / fulfills linear responsein the thermodynamic limit, meaning that the effect ofthe spin environment on the qudit is captured by thefirst two cumulants of the influence functional [36–38].For our model, the first cumulant vanishes and thus the reduced qudit dynamics is determined by the second cu-mulant, given by the force-force autocorrelation functionof the spin chain. This also means that one can constructan effective harmonic bath whose correlation function isthe same as that of the spin chain resulting in the samereduced qudit dynamics [37]. For our model, the effec-tive harmonic bath is characterized by a spectral densitywhich depends in general on the initial state, the ran-dom local fields and the spin-spin coupling g . For thespecific initial state considered here, the spectral densityof the effective harmonic bath is equal to the probabilitydistribution of twice the random local fields, and thus isindependent of g .Having established scaling collapse of the qudit vari-ance, we now turn our attention to dynamics of the spinchain, starting with the spin glass order parameter q EA ( t ) ≡ L − (cid:88) i (cid:104) ψ | σ z i ( t ) σ z i (0) | ψ (cid:105) . (5)Unlike the qudit variance, the spin glass order parameterdisplays marked drifts with system size (see insets of Fig.3(b,e)). The tendency of q EA ( t ) → γ , since γ controls the strength of a localtransverse field and thus governs the rate and magnitudeof a single spin’s precession. On reachable timescales t (cid:46) , the largest system size L = 96 has near perfectmemory of the initial state. This behavior is consistentwith our claim that the scaling of γ ∼ / √ L towardszero with increased system size is sufficiently fast thatthe system will flow to MBL for arbitrary γ , althoughproving MBL would require evolution to much later timesthan we can access.Though the usefulness of the influence functional ap-proach is restricted to the qudit, we should – by virtueof the fact that the initial spin dynamics are driven byinteractions with the qudit (for initial product states likethe super-Neel state we have chosen) – find that the spinobservables are linked to the qudit’s. The spin observ-ables should therefore enjoy a similar limiting behavioras γ ∝ L − / →
0. We indeed show this to be thecase within first order perturbation theory. In [24], weperform time-dependent perturbation theory using themethod of multiple scales. We solve for the time evolu-tion operator perturbatively by introducing new “inde-pendent” timescales t , t (cid:48) ≡ γt , t (cid:48)(cid:48) ≡ γ t , . . . , which allowfor control over secular terms growing with t . In thethermodynamic limit with scaled coupling, we find thatthe dynamics of the qudit are described perturbatively tofirst order up to time O (1 /γ ) (dotted lines in Fig. 3), pro-viding a complementary approach to the linear responsesolution from the influence functional formalism. Theperturbative calculation also demonstrates that spin ob-servables should exhibit similar gradual convergence toa single limit up to timescales t ∼ O (1 /γ ). Remark-ably, the connection between qudit variance and the spinglass order parameter is even more precise in this limit;they collapse to a single, universal curve in the thermody-namic limit upon scaling as ∆ /γ and (1 − q EA ) L/ (2 γ ),as seen in Fig. 3(a,b,d,e). Physically, this comes fromthe fact that a single perturbative excitation of the quditthrough the ˆ τ + + ˆ τ − component of H gives a single spinflip excitation of the spin chain through σ x j .Finally, we consider the entanglement entropy S A = − Tr [ ρ A log ρ A ] (6)between a contiguous half of the spins with the rest ofthe system, which is a defining feature in many bodylocalization. Here ρ A is the reduced density matrix ofhalf of the spin system, e.g., sites 1 through L/
2. As withthe previous two quantities, there appears to be a gradualconvergence of S A to a universal curve with increasing L ,although unlike the other observables, the entanglementdepends on the strength of the coupling prefactor γ .By turning off the Ising interaction g (dashed lines inFig. 3), we see that the dynamics of entanglement atshort times (cid:46) O (1) are unchanged – as predicted fromtime-dependent perturbation theory – while growth ofentanglement at intermediate times is dependent on this σ z i σ z i +1 interaction.These observations about the short-time dynamicshold for both weak and intermediate coupling, as seen inFig. 3. However, we can identify a slower timescale be-yond t (cid:46) O (1 /γ ) from first order perturbation theory, onwhich the Ising interactions start to play a role. In Fig.4, the same data is plotted upon rescaling the time by t/γ − . The observables are seen to roughly collapse forboth the weak and intermediate couplings and, for inter-mediate couplings, entanglement in particular shows in-teresting intermediate time behavior. While the collapse is imperfect, we note a few salient features. First, deepin the localized (weak coupling) regime, the spread of thequdit occupation, the growth of bipartite entanglemententropy, and the decay of the spin-glass order parameterappear to be arrested at long times. It is unclear whetherthe observables will continue to grow at later times, butour data leaves open the possibility that they saturateand that the asymptotic value may be system-size inde-pendent under the chosen scaling. Second, the dynamicsof the qudit appear to be correlated with dynamics ofthe spins, albeit with a slight time delay. Finally, in theintermediate coupling regime, the entanglement entropycontinues to grow at late times. For L = 16, there ap-pears to be a logarithmic growth over three decades inrescaled time (see Fig. 4g). The same may be true forthe L ≥
24, but we have insufficient data to decisivelyprove slow growth over several decades. As seen in Fig.4h, the period of potentially logarithmic growth coexistswith the period of finite occupation in the edge of thequdit spectrum (states n = 1 and d ), for which the high-frequency expansion yields all-to-all interactions (see Eq.2).It is unclear what drives the logarithmic behavior.When focusing on the bipartite entanglement entropy,two generic mechanisms have been studied in recentyears: the slow dephasing from a quench due to inter-actions between exponentially localized (quasi-local) op-erators [39, 40], and the linearly diverging semiclassicaltrajectories of the collective spin state [41] in long rangedinteracting spin systems. In the former case, it has beenfound that the slope of the logarithmic growth is indepen-dent of the strength of interactions [42]. This does notappear to be the case in our numerics, with the larger sys-tem sizes L ≥
24 ostensibly displaying log growth witha larger prefactor than in the L = 16 case. Moreover,there does not appear to be any logarithmic trend whenthe system is deep in the localized phase (see top row ofFig. 4). If conserved quasi-local operators do exist in thissystem, then our results would suggest that their local-ization lengths are strongly dependent on the coupling γ . Another possibility for the appearance of logarithmicgrowth of S A could come from the mediated all-to-allinteractions predicted in the effective Hamiltonian (Eq.2). In our qudit system, long ranged interactions beginto play a significant role when the extremal states of thequdit are occupied (see discussion in Sec. II). It was ar-gued that these mediated interactions are responsible forthe localization-delocalization transition upon decreasing d/ √ L , shown in Fig. 1. Consistent with this, we see sig-nificantly greater occupation in the extremal qudit statesfor intermediate couplings – where logarithmic growth isseen – compared to weak couplings (see Fig. 4(d,h)). Itis also clear that the slow growth of ∆ Q for intermedi-ate couplings is due in part to the slow growth in theoccupations of the | (cid:105) and | (cid:105) states.Regardless of the origin of slow growth, finite occu-pation at the extremes of the qudit spectrum implies a Ω t/ (2 π ) ∆ Q (a) L=8L=12L=16L=20 L=24L=48L=80 Ω t/ (2 π ) S A (b) FIG. 5. Dynamics for strong coupling, γ = 0 . / (cid:112) L/ (a) Qudit variance ∆ Q and (b) bipartite entanglement entropy S A . Results from ML-MCTDH are not included for ∆ Q asthey are not converged. For system sizes where the dynamicscan be computed exactly (dot-dashed lines), S A saturates thePage bound S A ∼ L/
2. The curves from ML-MCTDH (solidlines), corresponding to L ≥
24, saturate the bound set bythe number of single-particle functions in the second layer,log χ (dotted lines). departure from the Floquet regime. Our finite time nu-merics are unable to resolve whether this implies delo-calization. Should this mechanism give rise to a sharplocalization transition, it would possibly be of a differentcharacter from the extensively studied MBL transitionbased on ergodic grains thermalizing nearby insulatingregions through short range interactions [43–46]. B. Strong coupling (region II)
In the strong coupling regime (Fig. 5), our phase di-agram suggests that the system lies deep within thethermalizing phase for our available system sizes dueto strong delocalizing interactions between the spins in-duced by the central qudit. Our data is consistent withthis expectation, but we note two effects. First, de-spite the fact that spins thermalize, we observe that theasymptotic distribution of qudit occupations is nonuni-form, similar to the athermal qudit regime found in [12].In section D of [24] we introduce a phenomenological pic-ture to explain this based on random matrix theory, sug-gesting that it is rare for the qudit to make transitionsbetween widely separated states.Second, we note that due to the ergodic character ofthe dynamics in the strong-coupling thermalizing regime,the accurate treatment of the dynamics represents a sig-nificant challenge for the ML-MCTDH approach and can-not be converged for longer times [47]. This well knownlimitation of the ML-MCTDH method and other ten-sor network approaches is due to the following reason.Within the ML-MCTDH approach, the wave function ofthe system is represented in each layer by sums of Hartreeproducts, the total number of which is determined by thenumber N n of SPFs employed in a given layer n for each degree of freedom. For example, in the binary tree de-picted in Fig. 2, N SPFs are used in the second layerto represent each of the two parts of the spin chain re-sulting in ( N ) Hartree products that represent the spinsystem in the second layer. As a consequence, the en-tanglement entropy between the different constituents ofthe system is bounded by log N . However, for ergodicsystems, the entanglement entropy is extensive, and thus,starting from an uncorrelated state, the entanglement en-tropy grows and eventually exceeds the limit of log N .This implies that, for longer times, the wave function ofthe system cannot be represented accurately. The ap-plication of the ML-MCTDH formalism in the ergodicphase is thus restricted to short times. Therefore, the re-sults for the qudit variance depicted in Fig. 5 have beenobtained by exact diagonalization and the kernel polyno-mial method.Despite this limitation of ML-MCTDH, we are stillable to find signatures of thermalization by examiningthe dependence of the dynamics on N . In the rightpanel of Fig. 5, we see that the bipartite entanglemententropy is upper bounded by log N , corresponding toa maximal entropy state within our variational ansatz.We observe similarly strong dependence of q EA when in-creasing N , which drifts towards zero to indicate para-magnetic behavior in the spin chain. Other observables,such as the populations of the qudit levels cannot beconverged, implying that information about the ergodicstate is present, but limited. V. RESULTS FOR UNSCALED COUPLING
The dynamics of our system with scaled coupling, γ ∼ / √ L , is perhaps most interesting because it givesa well-defined thermodynamic limit for the qudit. How-ever, it is also important to understand the dynamicswhen the coupling is held fixed instead of being scaledby system size, corresponding to the dashed horizontalline in Fig. 1. Fixing the coupling strength may be easierto implement experimentally, for example in cavity QEDwhere the coupling is governed by the position-dependentelectric field strength. Doing so, however, means that wecan no longer easily separate dynamics occurring on dif-ferent timescales as in the previous section. Furthermore,we predict that in the thermodynamic limit this will even-tually result in thermalization, as the long-range interac-tions induced by the central qudit will eventually domi-nate at large enough times and system sizes. We choosethe fixed value γ = 0 . − Ω t/ (2 π ) . . . . . . − q E A (a) L=8L=12L=16L=24L=48L=96 − Ω t/ (2 π ) . . . . S A / L (b) − . . M I / L FIG. 6. Dynamics with unscaled coupling, γ = 0 . q EA (a) does not show significant system size dependence. How-ever, the bipartite entanglement entropy (b) , in addition tobeing subextensive – S A ∝ L α , 0 < α < L (cid:38) (inset) Mutual information MI ≡ I ( A, B )between two contiguous halves A and B of the spin chain (seemain text). The dynamics with fixed coupling, shown in Fig. 6,looks similar to the data at intermediate scaled couplingbut without as clear a separation of time scales or datacollapse. We note that it is harder to detect the sort oflogarithmically slow delocalization as seen in q EA in Fig.4f (cf. Fig. 6a). The story is the same with the quditvariance, which is similar to 1 − q EA when divided by L . However, with the fixed coupling, the entanglemententropy reflects a qualitative change in behavior at largeenough system sizes.We see that at times Ω2 π t ∼ × in the localizedphase (Fig. 6), the smaller system sizes L ≤
16 establisha subextensive amount of entanglement entropy. Numer-ics from ML-MCTDH seem to counter this trend, with S A continuing to grow slowly beyond this timescale. Therate of this growth increases with L , which is consistentwith it arising from stronger effective all-to-all interac-tions described by the high-frequency expansion (Eq. 2).It also appears to show strong system size dependenceat short times, where a subextensive amount of entan-glement is established. This should be contrasted withmodels of MBL without central coupling, in which theshort time behavior is system size independent.As a final note, we point out that the entanglemententropy at fixed γ is subextensive, such that S A /L ap-pears to be trending towards zero with increasing systemsize. This should be contrasted with mutual informationbetween the two halves of the spin chain, I ( A, B ) = S ( A ) + S ( B ) − S ( A ∪ B ) , which is extensive. In Ref. [12], mutual information wasused as a proxy for entanglement between the two halves of the spin chain, as it nominally removes “unimportant”entanglement with the central qudit. However, since en-tanglement must be subextensive – indeed, system sizeindependent – in the MBL phase, our data indicate thatentanglement entropy is a better metric than mutual in-formation for capturing this. Our initial expectation wasthat mutual information would become subextensive atlarger system size, but the results obtained with the ML-MCTDH method rule out that possibility. VI. CONCLUSIONS
In this paper we have studied the dynamical behaviorof a qudit coupled to a disordered, interacting bath ofup to L = 96 spins-1 /
2, which altogether can exhibitlocalization at strong disorder. Using a combination ofexact propagation methods and the tensor network-basedML-MCTDH approach, we find evidence of qualitativelydifferent dynamical signatures in local observables suchas the spin glass order of the bath and the qudit variance,consistent with a rough phase diagram (Fig. 1). Mostnotably, we find hints of logarithmically slow decay oflocalization near the onset of all-to-all interactions in thebath. This behavior was found to occur after timescales t ∼ O (1 /γ ) where γ is the qudit-spin bath coupling.The behavior of the qudit observed here is, we believe,not specific to this model. Our conclusions should ap-ply equally well to the cases of a cavity photon withrescaled raising/lowering operators a † → ( N ) − / a † orcentral spin- S systems with operators rescaled as ˆ S → ( S ( S − − / ˆ S . The feature of these systems is thatthe fundamental commutation relation between the rais-ing and lowering operators vanishes in the limit of large S or large N . This fact allows for exact cancellationbetween processes that raise or lower the qudit state.However, this mechanism only serves to protect localiza-tion for sufficiently large “magnetic field” Ω; it is unclearhow these systems interpolate between the Ω = 0 limitand the Ω > | g | , | h i | , . . . limit. We note additionally thatthe limitations of ML-MCTDH for these types of cen-trally coupled systems with many-body interacting bathsin the strong coupling regime requires more clarification.Such clarifications may be necessary to extend the ef-fectiveness of the method into the thermalizing regimeon the left side of the phase diagram 1, which remainsnumerically inaccessible and thus poorly understood. Acknowledgements
This work was performed with support from the Na-tional Science Foundation through award number DMR-1945529 (MHK), the Welch Foundation through awardnumber AT-2036-20200401 (MHK), and the German Re-search Foundation (DFG) through IRTG 2079. This re-search used resources of the National Energy ResearchScientific Computing Center, a U.S. Department of En-ergy Office of Science User Facility operated under Con-tract No. DE-AC02-05CH11231. Furthermore, support0by the state of Baden-W¨urttemberg through bwHPC and the DFG through grant no. INST 40/575-1 FUGG (JUS-TUS 2 cluster) is gratefully acknowledged. [1] A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, andR. J. Schoelkopf, Phys. Rev. A , 062320 (2004).[2] E. J. Davis, G. Bentsen, L. Homeier, T. Li, and M. H.Schleier-Smith, Phys. Rev. Lett. , 010405 (2019).[3] K. Xu, J.-J. Chen, Y. Zeng, Y.-R. Zhang, C. Song,W. Liu, Q. Guo, P. Zhang, D. Xu, H. Deng, K. Huang,H. Wang, X. Zhu, D. Zheng, and H. Fan, Phys. Rev.Lett. , 050507 (2018).[4] J. Z. Imbrie, Journal of Statistical Physics , 998(2016).[5] D. A. Huse, R. Nandkishore, V. Oganesyan, A. Pal, andS. L. Sondhi, Phys. Rev. B , 014206 (2013).[6] V. Khemani, Ph.D. thesis (2016).[7] N. Y. Yao, C. R. Laumann, S. Gopalakrishnan, M. Knap,M. M¨uller, E. A. Demler, and M. D. Lukin, Phys. Rev.Lett. , 243002 (2014).[8] A. O. Maksymov and A. L. Burin, Phys. Rev. B ,024201 (2020).[9] R. Nandkishore, S. Gopalakrishnan, and D. A. Huse,Phys. Rev. B , 064203 (2014).[10] D. A. Abanin, W. D. Roeck, and F. Huveneers, Annalsof Physics , 1 (2016).[11] P. Ponte, Z. Papi´c, F. Huveneers, and D. A. Abanin,Phys. Rev. Lett. , 140401 (2015).[12] N. Ng and M. Kolodrubetz, Phys. Rev. Lett. , 240402(2019).[13] P. Ponte, C. R. Laumann, D. A. Huse, andA. Chandran, Phil. Trans. R. Soc. A. (2017),10.1098/rsta.2016.0428.[14] D. Hetterich, N. Y. Yao, M. Serbyn, F. Pollmann, andB. Trauzettel, Phys. Rev. B , 161122 (2018).[15] H. Wang and M. Thoss, J. Chem. Phys. , 1289 (2003),https://doi.org/10.1063/1.1580111.[16] U. Manthe, J. Chem. Phys. , 164116 (2008),https://doi.org/10.1063/1.2902982.[17] H.-D. Meyer, F. Gatti, and G. Worth,Multidimensional Quantum Dynamics (Wiley-VCH,Weinheim, 2009).[18] O. Vendrell and H.-D. Meyer, J. Chem. Phys. ,044135 (2011), https://doi.org/10.1063/1.3535541.[19] H. Wang, The Journal of Physical Chem-istry A , 7951 (2015), pMID: 26020459,https://doi.org/10.1021/acs.jpca.5b03256.[20] M. Schreiber, S. S. Hodgman, P. Bordia, H. P.L¨uschen, M. H. Fischer, R. Vosk, E. Altman,U. Schneider, and I. Bloch, Science , 842 (2015),https://science.sciencemag.org/content/349/6250/842.full.pdf.[21] P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas,A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro,A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jef-frey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quin-tana, D. Sank, A. Vainsencher, J. Wenner, T. White,H. Neven, D. G. Angelakis, and J. Martinis, Science , 1175 (2017).[22] P. Sierant, K. Biedro´n, G. Morigi, and J. Zakrzewski,SciPost Phys. , 8 (2019).[23] This estimate is subject to strong finite size effects. [24] See Supplemental Material at URL for discussion on thecumulant expansion, an ansatz for the steady state ofthe thermalizing phase, and perturbation theory usingmultiple scales.[25] S. Paeckel, T. K¨ohler, A. Swoboda, S. R. Manmana,U. Schollw¨ock, and C. Hubig, Annals of Physics ,167998 (2019).[26] A. Lubatsch and R. Frank, The European Physical Jour-nal B , 1 (2019).[27] G. Worth, M. H. Beck, A. J¨ackle, O. Vendrell, and H.-D.Meyer, The MCTDH Package, Version 8.5.11 (2019), seehttp://mctdh.uni-hd.de.[28] H. Wang and J. Shao, J. Chem. Phys. , 22A504(2012), https://doi.org/10.1063/1.4732808.[29] H.-D. Meyer, U. Manthe, and L. Cederbaum, ChemicalPhysics Letters , 73 (1990).[30] U. Manthe, H.-D. Meyer, and L. S. Cederbaum, J. Comp.Phys. , 3199 (1992).[31] M. Beck, A. J¨ackle, G. Worth, and H.-D. Meyer, PhysicsReports , 1 (2000).[32] H.-D. Meyer, Wiley Interdisciplinary Reviews: Compu-tational Molecular Science , 351 (2012).[33] H. Wang and M. Thoss, J. Chem. Phys. , 024114(2009), https://doi.org/10.1063/1.3173823.[34] U. Weiss, Quantum Dissipative Systems (World Scien-tific, 2012).[35] M. Kac, G. E. Uhlenbeck, and P. C. Hem-mer, Journal of Mathematical Physics , 216 (1963),https://doi.org/10.1063/1.1703946.[36] R. Feynman and F. Vernon, Annals of Physics , 118(1963).[37] N. Makri, The Journal of Physical Chemistry B ,2823 (1999), https://doi.org/10.1021/jp9847540.[38] H. Wang and M. Thoss, J. Phys. Chem. A , 10369(2007).[39] M. ˇZnidariˇc, T. c. v. Prosen, and P. Prelovˇsek, Phys.Rev. B , 064426 (2008).[40] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.Mod. Phys. , 021001 (2019).[41] A. Lerose and S. Pappalardi, Phys. Rev. Research ,012041 (2020).[42] J. A. Kj¨all, J. H. Bardarson, and F. Pollmann, Phys.Rev. Lett. , 107204 (2014).[43] R. Vosk and E. Altman, Phys. Rev. Lett. , 067204(2013).[44] P. T. Dumitrescu, R. Vasseur, and A. C. Potter, Phys.Rev. Lett. , 110604 (2017).[45] A. Goremykina, R. Vasseur, and M. Serbyn, Phys. Rev.Lett. , 040601 (2019).[46] A. Morningstar, D. A. Huse, and J. Z. Imbrie, Phys.Rev. B , 125134 (2020).[47] T. Westermann and U. Manthe, J. Chem. Phys. ,204116 (2012), https://doi.org/10.1063/1.4720567.[48] J. Shao and P. H¨anggi, Phys. Rev. Lett. , 5710 (1998).[49] H. Wang and M. Thoss, J. Phys. Chem. A , 10369(2007).[50] A. Lucas, “Quantum many-body dynamics on the star graph,” (2019), arXiv:1903.01468 [cond-mat.str-el].[51] P. Hauke and M. Heyl, Phys. Rev. B , 134204 (2015). Supplemental Material for “Localization dynamics in a centrally coupled system”
Appendix A: Linear response property in the thermodynamics limit
In the limit L → ∞ the spin chain can be seen as a macroscopic environment allowing for a different viewpointon the dynamics of the qudit. Environments consisting of independent (i.e. noninteracting) degrees of freedom andappropriate scaling of the coupling γ fulfil linear response theory, i.e. the influence of the environment on the qudit isfully characterized by the force autocorrelation function of the bath allowing for the description of the influence of theenvironment in terms of a bath of harmonic oscillators with an effective spectral density.[37] This idea was used, forexample, to describe the dynamics of a two-level system coupled to a spin bath consisting of independent spins[37, 48]and to a bath of anharmonic vibrational degrees of freedom[49]. For non-interacting environments this was provenby showing that all but the leading order term in the cumulant expansion of the influence functional vanish in thethermodynamic limit. Despite the fact that the spin chain considered in our work consists of interacting spins, onecan show that the linear response property also holds for this model, and thus, the influence of the spin chain on thequdit can be described by a bath of harmonic oscillators with an effective spectral density.To this end, we follow the derivation given in Ref. [37] and extend it to the interacting spin chain considered here.We assume that initially there are no correlations between different spins and the spin chain is in an eigenstate of all { σ zi } . In order to evaluate the different terms in the cumulant expansion[37] the time dependent operatorˆ f ( t ) = − γ L (cid:88) i =1 e iH t σ xi e − iH t (cid:124) (cid:123)(cid:122) (cid:125) σ xi ( t ) (A1)is required, which corresponds to the bath part of the system-bath interaction and describes the force exerted on thesystem due to its interaction with the environment. Here, H is the Hamiltonian of the isolated spin chain. The timeevolved operators σ xi ( t ) = σ + i ( t ) + σ − i ( t ) can be calculated analytically yielding σ + i ( t ) = e ihξ i t (cid:20) [ − σ zi − σ zi +1 ] + [ + σ zi − σ zi +1 ] cos(4 gt ) + i [ σ zi − + σ zi +1 ] sin(4 gt ) (cid:21) σ + i = e ihξ i t ˆ ϕ i ( t ) σ + i , (A2)and σ − i ( t ) = ( σ + i ( t )) † . Here σ + i and σ − i are the spin raising and lowering operators, respectively. Using this equation,one can show that (cid:2) σ xi ( t ) , σ xj ( t (cid:48) ) (cid:3) = 0 ∀ i, j with | i − j | ≥ , (A3)This follows from the fact that ˆ ϕ i ( t ) only involves spin operators on sites i − i + 1.All terms in the cumulant expansion of the influence functional can be expressed[37] in terms of N -time correlationfunctions defined as C ( N ) ( t , ..., t N ) = (cid:104) ˆ f ( t ) ... ˆ f ( t N ) (cid:105) , (A4)where (cid:104) ... (cid:105) denotes the expectation value with respect to the initial state of the environment. Initially, the spin chainis in an eigenstate of all { σ zi } , and thus, one finds that C (1) = 0 , (A5)Consequently, the first term, and by extension all odd order terms in the expansion, vanishes. The second order termcan be expressed in terms of the two-point correlation function C (2) ( t , t ) = γ (cid:88) i (cid:88) j (cid:104) σ xi ( t ) σ xj ( t ) (cid:105) . (A6)It is straightforward to check that the expectation value is zero for all i (cid:54) = j , and thus, the two-time correlationfunction reduces to C (2) ( t , t ) = γ L (cid:88) i =1 (cid:104) σ xi ( t ) σ xi ( t ) (cid:105) . (A7)3Since (cid:104) σ xi ( t ) σ xi ( t ) (cid:105) is finite but not zero in general, the sum diverges in the thermodynamic limit L → ∞ unless γ ∼ / √ L , which we will assume in the following by setting γ = γ (cid:112) L / L . The next non-vanishing term in thecumulant expansion is the fourth term which involves the fourth order correlation function C (4) ( t , t , t , t ) = γ (cid:88) i,j,k,l (cid:104) σ xi ( t ) σ xj ( t ) σ xk ( t ) σ xl ( t ) (cid:105) . (A8)Since ˆ ϕ i ( t ) does not change the initial state of the spin chain, there can only be up to two different indices in theexpectation value. Thus, the four-time correlation function can be written as C (4) ( t , t , t , t ) = γ (cid:88) i,j (cid:104) σ xi ( t ) σ xi ( t ) σ xj ( t ) σ xj ( t ) (cid:105) + γ (cid:88) i,ji (cid:54) = j (cid:104) σ xi ( t ) σ xj ( t ) σ xi ( t ) σ xj ( t ) (cid:105) + γ (cid:88) i,ji (cid:54) = j (cid:104) σ xi ( t ) σ xj ( t ) σ xj ( t ) σ xi ( t ) (cid:105) . (A9)In the following we will discuss only the first term on the right hand side. The other two terms can be treatedequivalently. The double sum can be decomposed as γ (cid:88) i,j (cid:104) σ xi ( t ) σ xi ( t ) σ xj ( t ) σ xj ( t ) (cid:105) = γ (cid:88) i (cid:104) σ xi ( t ) σ xi ( t ) σ xi ( t ) σ xi ( t ) (cid:105) + γ (cid:88) i (cid:104) σ xi ( t ) σ xi ( t ) σ xi +1 ( t ) σ xi +1 ( t ) (cid:105) + γ (cid:88) i (cid:104) σ xi +1 ( t ) σ xi +1 ( t ) σ xi ( t ) σ xi ( t ) (cid:105) + γ L (cid:88) i,j | i − j |≥ (cid:104) σ xi ( t ) σ xi ( t ) σ xj ( t ) σ xj ( t ) (cid:105) . (A10)The expectation values in the first three sums on the right hand side are bounded by a constant. If the coupling γ isscaled as γ = γ (cid:112) L / L , these terms vanish as L → ∞ , and thus, we conclude that γ L (cid:88) i,j (cid:104) σ xi ( t ) σ xi ( t ) σ xj ( t ) σ xj ( t ) (cid:105) = γ L (cid:88) i,j | i − j |≥ (cid:104) σ xi ( t ) σ xi ( t ) σ xj ( t ) σ xj ( t ) (cid:105) + O (cid:18) L (cid:19) . (A11)Because the operators in the expectation value act on different Hilbert spaces and the initial state factorizes theexpectation values can be factorized as γ L (cid:88) i,j | i − j |≥ (cid:104) σ xi ( t ) σ xi ( t ) σ xj ( t ) σ xj ( t ) (cid:105) = γ L (cid:88) i,j | i − j |≥ (cid:104) σ xi ( t ) σ xi ( t ) (cid:105) (cid:104) σ xj ( t ) σ xj ( t ) (cid:105) (A12)Adding the terms for i = j and | i − j | = 1 to the double sum on the right hand side gives an error of O ( / L ), andthus, one can write γ (cid:88) i,j (cid:104) σ xi ( t ) σ xi ( t ) σ xj ( t ) σ xj ( t ) (cid:105) = γ (cid:88) i,j (cid:104) σ xi ( t ) σ xi ( t ) (cid:105) (cid:104) σ xj ( t ) σ xj ( t ) (cid:105) + O (cid:18) L (cid:19) (A13)= C (2) ( t , t ) C (2) ( t , t ) + O (cid:18) L (cid:19) , (A14)where we have identified the two-time correlation functions. With this, we finally conclude thatlim L →∞ C (4) ( t , t , t , t ) = lim L →∞ C (2) ( t , t ) C (2) ( t , t )+ C (2) ( t , t ) C (2) ( t , t )+ C (2) ( t , t ) C (2) ( t , t ) . (A15)4Using this one finds that the fourth order term in the cumulant expansion in [37] vanishes in the thermodynamic limit.In a similar way one can show that all higher order terms in the expansion vanish, proving that in the thermodynamiclimit the influence functional is completely characterized by the force autocorrelation function. Thus, one can constructa bath of harmonic oscillators with an effective spectral density resulting in the same influence functional.As a last step we show that for the choice of the ’super-Neel’ state | E SN (cid:105) = |↑↑↓↓ ... (cid:105) as initial state the forceautocorrelation function does not depend on the spin-spin interaction g . The force autocorrelation function is definedas (cid:104) ˆ f ( t ) ˆ f ( t ) (cid:105) = γ L (cid:88) i (cid:104) E SN | σ xi ( t ) σ xi ( t ) | E SN (cid:105) , (A16)where the time-dependent operators are given in equation (A2). For the ”super-Neel” state it follows that σ zi − σ zi +1 | E SN (cid:105) = − | E SN (cid:105) , (cid:0) σ zi − + σ zi +1 (cid:1) | E SN (cid:105) = 0 , (A17)holds for all i since the spins at site i − i + 1 are always antiparallel. Thus, the action of ˆ ϕ i ( t ) is independentof the index i and gives ˆ ϕ i ( t ) | E SN (cid:105) = | E SN (cid:105) . (A18)Since the spin-spin interaction enters only via ˆ ϕ i ( t ), we find that the force autocorrelation function, and conse-quently the parameters of the effective bath harmonic oscillators, are independent of the spin-spin interaction g . Thecorresponding effective spectral density can be calculated from the force autocorrelation[37] yielding J eff ( ω ) = π h χ [ − h, h ] ( ω ) , (A19)where χ I is the characteristic function of the interval I , i.e. χ I ( ω ) = 1 if ω ∈ I and 0 else. Thus, in the thermodynamicslimit a bath of harmonic oscillators with this spectral density gives rise to the same qudit dynamics as the spin chainin the ”Super-Neel” state. Appendix B: Structure of wavefunction
In addition to modifying the entanglement dynamics at short times, the star-like geometry of this system (depictedin the inset of Figure 2) should render the concept of locality meaningless. Indeed, from the point of view of operatordynamics, operators for the qudit should immediately spread to O ( L ) sites after O (1) time[50]. Instead we proposeto analyze the structure of eigenstates in the Hilbert space of the uncoupled ( γ = 0) Hamiltonian. In the spirit ofthe current understanding of MBL, where eigenstates are weak deformations of the unperturbed system, we parsewavefunctions in the product basis | s (cid:105) of spins and qudit states – | z (cid:105) ⊗ | n (cid:105) with z i ∈ {↑ , ↓} , ∀ i = 1 , . . . , L and n = 1 , . . . , d . We quantify “deformations” through a notion of Hilbert space distance as D ( | z , n (cid:105) , | z (cid:48) , n (cid:48) (cid:105) ) = max ( D H ( z , z (cid:48) ) , | n − n (cid:48) | ) , (B1)with D H being the Hamming distance between the bitstrings z and z (cid:48) (here, up-spins are “1” and down-spins are“0”). Intuitively, this measures the minimum number of times the transverse perturbation must be applied to connecttwo states in the Hilbert space. For eigenstates, the origin ( ≡ | ψ ref (cid:105) ) is taken to be the product state with the largestweight while for time evolution, the origin is taken to be the initial state before the quench. Other product statescan then be grouped according to their distance D from | ψ ref (cid:105) . For each D , we calculate the distribution of expansioncoefficients | (cid:104) s | ψ (cid:105) | over disorder realizations and equidistant product states {| s (cid:105) | D ( | ψ (cid:105) , | s (cid:105) ) = x } . In the nonergodicphase, these coefficients are suppressed by a factor of ( γ/g ) as D increases (Figure 7(a,b)). One can then regard thewavefunction as being exponentially localized in Hilbert space. This can be observed also in non-centrally coupledmodels of MBL, such as the disordered Ising chain with next-nearest neighbor interactions [42]. In contrast, there isno such Hilbert space localization at large enough γ in the ergodic phase (inset of Figure 5b). This is corroboratedby the average Hilbert space distance (cid:104)D ( t ) (cid:105) , as measured from | ψ ref (cid:105) . This quantity has been noted by Hauke andHeyl [51] to saturate to L/ (cid:104)D ( t ) (cid:105) (noted by [51]) approximately holds (see top panel ofFigure 7c), i.e. (cid:104)D ( t ) (cid:105) ≥ L (1 − q EA ) / D we have chosen.5 (a) - - - - L =
12, d = γ = ϵ = 〈 s | ψ〉 Larger log ( γ / g ) = (b) - - - - L =
12, d = γ = = log 〈 s | ψ〉 = - - - - - γ = ϵ = (c) t / ℏ 〈 D ( t ) 〉 / L 〈 D 〉 / L - ( - q E A ) / × - γ = γ = γ = γ = = = FIG. 7. Distributions of wavefunction coefficients in the unperturbed basis {| s (cid:105)} , parsed by the Hilbert space distance D separating | s (cid:105) and a reference state | s ref (cid:105) . The colors of the curves lighten for increasing D , while the dashed vertical lines areguides to the eye for where the expected peak locations should the coefficients decay as ∝ ( γ/g ) D . (a) Distributions in the teneigenstates closest to the middle of the energy spectrum, choosing | s ref (cid:105) to be the unperturbed state with largest weight. (b) Distribution of coefficients for the system at t = 10 , evolving from | ψ (0) (cid:105) = | s ref (cid:105)| ( d + 1) / (cid:105) with | s ref (cid:105) being the super-Neelstate. (c) Average Hilbert space distance (cid:104) D ( t ) (cid:105) from the super-Neel state (bottom) and in comparison with the spin glassorder parameter q EA (top), for γ = 0 . , . , . , . L = 8(dashed) , In the top panel of Figure 7c, the periods where D /L > (1 − q EA ) / ∼ /γ , versus the slower decay of q EA . The latter proceeds on slower timescalesthrough a combination of qudit-spin flip transitions and the Ising interaction gσ z i σ z i +1 . This interpretation of the spinglass order parameter bring new meaning to the results from ML-MCTDH. At least up to intermediate times, thelocalization length of the wavefunction in Hilbert space is stable up to L = 96 when γ ≈ . Appendix C: Trivial limit and convergence of qudit variance
The quench setup we examine in this paper – in which we prepare the full system in an eigenstate of the γ = 0Hamiltonian – results in certain behaviors in the γ → only generator ofdynamics in both the spins and the qudit at short times. We demonstrate this limit by plotting observables rescaledby γ − in Figure 8. This scaling at small γ converges the dynamics at short times up to t ≤ . At small enoughcoupling, a plateau begins to appear at t ∼ . We see hints of this in the MCTDH data for mutual information(Figure 4c), for example, with the establishment of a plateau for L = 48 and 80 at times 10 ≤ t (cid:46) . This behaviorcontrasts with the immediate increase of MI for smaller system sizes at t ∼ (a) - γ - ( - q ) (b) - γ - Δ Q γ = γ = γ = γ = γ = (c) - γ - M I FIG. 8. Dynamics from the “super-Neel” state for L = 8 across different γ in the localized phase. Rescaling the (a) spin glassorder parameter, (b) qudit variance, and (c) mutual information by the coupling collapses the dynamics at short to intermediatetimes. The asymptotic behavior in each case should be interpreted as defining a trivial limit analogous to Anderson localizationvis-a-vis MBL. We take advantage of the vanishing coupling, which is scaled to zero with system size, to systematically re-6construct the dynamics using the method of multiple scales. To that end, we solve for the time evolution op-erator, artificially introducing new ‘independent’ timescales t , t (cid:48) ≡ γt , t (cid:48)(cid:48) ≡ γ t , . . . , which allow for con-trol over secular terms growing unboundedly as ∼ t . Formally, the time development operator is expanded as U ( t ) ≡ U ( t, t (cid:48) , t (cid:48)(cid:48) , . . . ) + γU ( t, t (cid:48) , t (cid:48)(cid:48) , . . . ) + γ U ( t, t (cid:48) , t (cid:48)(cid:48) , . . . ) + . . . and we shall solve for the full evolution order-by-order. We note that we have made an important assumption that the only timescales of interest are ∼ γ − n . However,we are not interested in describing the dynamics for all t at arbitrary γ and L , but only up to the t ∼ γ − as γ → ddt U ( t ) = − i ( H + Ω τ z + γH τ x ) U ( t ) , where we define τ z = d (cid:88) n =1 n | n (cid:105)(cid:104) n | , τ + = d − (cid:88) n =1 | n + 1 (cid:105)(cid:104) n | , and τ − = ( τ + ) † . From these we construct τ x = τ + + τ − and τ y = − iτ + + iτ − . With the addition of the new timescales, the time derivative now becomes ddt U = (cid:18) ∂∂t U (cid:19) + γ (cid:18) ∂∂t (cid:48) U + ∂∂t U (cid:19) + γ (cid:18) ∂∂t (cid:48)(cid:48) U + ∂∂t (cid:48) U + ∂∂t U (cid:19) + . . . At the zeroth order, the equation of motion and its solution are ∂∂t U = − i ( H + Ω τ z ) U = ⇒ U = e − i ( H +Ω τ z ) t U int0 ( t (cid:48) , t (cid:48)(cid:48) , . . . ) such that U int0 (0 , , . . . ) = 1 . At first order, ∂∂t (cid:48) U + ∂∂t U = − i ( H + Ω τ z ) U − iH τ x U ∂∂t (cid:48) U int0 + ∂∂t U int1 = − ie i ( H +Ω τ z ) t H τ x e − i ( H +Ω τ z ) t V , (C1)where we let U = e − iH t U int1 . Note at this point that the first term on the LHS is independent of t . Its contribution to U int1 would be proportional to t and thus secular. Should there also exist secular terms on the RHS (i.e., independentof t ), U int0 should be chosen to offset it. Otherwise, it must be independent of t (cid:48) , i.e. U int0 ≡ V ( t (cid:48)(cid:48) , . . . ). If these secularterms do not exist at all orders, then the multiple scales result would be completely equivalent to the Dyson series inthe interaction picture. In our disordered system, we must be careful of secularity and near-secularity. The former, inwhich two states linked by the perturbation are exactly degenerate, occurs with zero probability since the local fieldmust have a value such that Ω = ± h i + g ( σ z i − + σ z i +1 )). More likely is the scenario of near-degeneracies, which atthis order can lead to U int1 growing arbitrarily large after an arbitrarily long time. Such terms make the expansion of U uncontrolled at long times. We can, however, absorb near-secular behavior into U int0 .Define A ( t (cid:48) ) = (cid:88) | a (cid:105) , | b (cid:105)| E a − E b | < (cid:104) a | H τ x | b (cid:105) exp (cid:18) i ∆ E ab γ t (cid:48) (cid:19) | a (cid:105)(cid:104) b | . To regulate the secular part of (C1), we must have ∂∂t (cid:48) U int0 = − iA ( t (cid:48) ) U int0 . This makes U int0 ( t (cid:48) , t (cid:48)(cid:48) , . . . ) = exp − i t (cid:48) (cid:90) A ( τ ) dτ V ( t (cid:48)(cid:48) , . . . ) . The unknown function V will be solved for at higher orders. The argument of the exponential should always becomplex, since A ( t (cid:48) ) is Hermitian. Thus these resonant terms will not cause U int0 to have unbounded norm, and theperturbative expansion for U remains valid. However, our ability to regulate the secularity in this way should not betaken as a statement on the dynamics being localized. Instead, it implies that more careful consideration of U int0 isnecessary to understand if resonances are able to cause delocalization. At present, all we need is to examine if we cansafely neglect exp( − i (cid:82) γt A ) if we scale γ → A as roughly being theadjacency matrix for states in Hilbert space, where two states are connected by an edge if they are resonant. It isknown that the eigenvalues of adjacency matrices are bounded above by the maximum number of edges connecting to7a vertex in the graph. A naive upper bound for our system is 2 L , which is the number of states that can be reachedby applying the perturbation on to the eigenstates of the unperturbed system. Should the spectrum of A saturatethis bound, the matrix exponential will contain time dependence going as exp( − i γtL ) and scaling γ ∝ L − / willnot remove the correction factor in U int0 . Regardless, this will not pose a problem in our model for the parametersand initial super-Neel state we have chosen.Having found the lowest order approximation for U int0 , we are now left with the nonsecular terms. These can bestraightforwardly used to solve for U int1 . Let ( H τ x ) reg be the regular version of the perturbation H τ x with resonantmatrix elements removed. The equation of motion becomes ∂∂t U int1 = − i e i ( H +Ω τ z ) t ( H τ x ) reg e − i ( H +Ω τ z ) t U int0 ( t (cid:48) , . . . )= ⇒ U = − i e − i ( H +Ω τ z ) t iV ( t (cid:48) , t (cid:48)(cid:48) , . . . ) + t (cid:90) e i ( H +Ω τ z ) τ ( H τ x ) reg e − i ( H +Ω τ z ) τ dτ U int0 ( t (cid:48) , . . . ) While we can in principle keep going to higher orders, we will stop here and discuss the dynamics with scaledqudit-spin chain coupling. The results we have obtained so far allow us to accurately describe time evolution up to atimescale t ∼ O (1 /γ ). Should we keep decreasing γ , then all the artificial times t (cid:48) , t (cid:48)(cid:48) , . . . , will tend to zero withoutaffecting the physical time t . Using the initial condition for the unknown functions can then give us closed formexpressions for U . For example, we approximate U ( t ) ≈ e − i ( H +Ω τ z ) t + γ e − i ( H +Ω τ z ) t − i t (cid:90) e i ( H +Ω τ z ) τ ( H τ x ) reg e − i ( H +Ω τ z ) τ dτ . This is essentially what one finds in usual perturbation theory, except we now have more knowledge of its convergenceproperties.We can gain some analytical understanding of the qudit variance and the spin glass order parameter from thisapproximation to the propagator.Numerics show that the dynamics of qudit variance ∆ Q = (cid:10) τ ( t ) (cid:11) − (cid:104) τ ( t ) (cid:105) comes mostly from the first term, asthe second term is essentially constant. We calculate ddt ∆ Q ≈ (cid:10)(cid:8) τ z , ddt τ z (cid:9)(cid:11) , where ddt τ z = γH τ y . ddt ∆ Q ≈ γ (cid:104) ψ | { τ z ( t ) , H ( t ) τ y ( t ) } | ψ (cid:105) = γ (cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:58) (cid:68)(cid:110) τ z (0) , H (0)1 τ y (0) (cid:111)(cid:69) + γ (cid:68)(cid:110) τ z (1) , H (0)1 τ y (0) (cid:111)(cid:69) + γ (cid:68)(cid:110) τ z (0) , (cid:16) H (1)1 τ y (0) + H (0)1 τ y (1) (cid:17)(cid:111)(cid:69) . The first term is zero since the operator is off-diagonal. The last term must also vanish since it is not invariant withrespect to redefinition of τ z , e.g. changing its spectrum from (0 , . . . , d −
1) to (1 , . . . , d ) by adding a constant termto the Hamiltonian. Indeed one can check that it is exactly cancelled by the − (cid:104) τ z (cid:105) term we have neglected in ourapproximation of the qudit variance. A calculation of the remaining term shows ddt ∆ Q ≈ γ (cid:88) i (cid:88) ± sin( t ∆ E ± i )∆ E ± i (cid:68) ψ ± i (cid:12)(cid:12)(cid:12) ( σ x i τ x ) reg (cid:12)(cid:12)(cid:12) ψ (cid:69) , (C2)consistent with usual perturbation theory. Upon disorder averaging, we see that∆ Q ≈ γ L t (cid:90) dτ (cid:88) ± sin( τ ∆ E ± i )∆ E ± i (cid:68) ψ ± i (cid:12)(cid:12)(cid:12) ( σ x i τ x ) reg (cid:12)(cid:12)(cid:12) ψ (cid:69) . Thus the qudit variance – along with other qudit observables such as the population – converge to a single curveupon scaling the coupling as γ ∝ / √ L . Convergence towards this expression should be expected up to time t ∼ O ( γ − ) ∝ O ( √ L ). For dynamics from the super-Neel state as we study here, the energy difference with spin i flippedis ∆ E i = ± h i . There are no resonances for our chosen values of h i ∈ [ − . , .
3] and Ω ≈ .
93, so we have exactly∆ Q ≈ γ (cid:88) i (cid:88) ± − cos ( t (2 h i s i ± Ω))(2 h i s i ± Ω) ddt ∆ Q ≈ γ L (cid:18) Si[ t (2 h + Ω)] + Si[ t (2 h − Ω)] h (cid:19) , t ) is the sine integral. For fixed γ at large enough L , this expression will violate the bound on the quditvariance, ( d − /
12, when all the states of the qudit are equally populated. Thus higher order terms are necessaryto prevent this unphysical outcome.We can similarly look at the spin glass order parameter, and find that for the super-Neel state, ddt (1 − q ) = 4 γ L (cid:88) i (cid:88) ± sin( t ∆ E ± i )∆ E ± i . This shows the surprising fact that the qudit variance (cf. (C2)) and spin glass order are linearly related to each otherin this limit. This motivates the scaling we take in plotting the results in Figures 3 and 4 in the main text.
Appendix D: Ansatz for qudit populations in the thermalized phase
We consider eigenstate thermalization in the sense that (cid:104) ψ | A | ψ (cid:105) = Z − Tr (cid:2) e − βH A (cid:3) , where β is the temperature reproducing the same energy (cid:104) ψ | H | ψ (cid:105) . When we time evolve from an initial state sitting atthe middle of the many-body spectrum, | ψ (0) (cid:105) = | ↑↑↓↓ . . . (cid:105) (cid:12)(cid:12) d − (cid:11) , we should consider infinite temperature averages,i.e. (cid:104) ψ ( t ) | A | ψ ( t ) (cid:105) = d − Tr qd (cid:0) − L Tr S A (cid:1) . Therefore by setting A = | n (cid:105)(cid:104) n | for n = (1 − d ) / , . . . , ( d − /
2, we should expect a uniform occupation over all d levels of the qudit. We do not observe this in the delocalized phase; instead, the occupations of the qudit seem tosaturate to the distributions shown in Fig. 9 (averaging over ∼
10 realizations of disorder). ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ◆ - - - n p n Qudit occupationsScaled coupling γ = / L / ● L = (
10 realizations ) ■ L = (
10 realizations ) ◆ L = ( ) ● ● ● ●■ ■ ■ ■ ◆ ◆ ◆ ◆ n p n Symmetrized qudit occupationsScaled coupling γ = / L / ● L = (
10 realizations ) ■ L = (
10 realizations ) ◆ L = ( ) σ = σ = σ = FIG. 9. Occupations of the qudit at t ∼ g − . On the right, the occupations have been vertically spaced by 0 .
05 for clarity.
On the left, we notice that the occupations appear to be asymmetric about the middle state. This is possibly due tothe small number of disorder averages, but may also be affected by finite size effects or the rather special super-Neelinitial state. As a proxy for the infinite sample limit, we symmetrize the occupations and note that a one parameterGaussian ansatz, p n ∝ exp (cid:0) − n / (2 σ ) (cid:1) , fits well.We propose a model to reproduce these observations based on the partial diagonalization achieved by invokingFloquet’s theorem. By going into the rotating frame through the transformation exp ( − i Ωˆ nt ), the Hamiltonianbecomes time-periodic and we can factorize the time evolution in the “lab frame” as U ( t ) = e − i Ωˆ nt e − iK rot ( t ) e − iH effrot t e iK rot (0) . We had previously proven that eigenstates can be written in the form | E (cid:105) = e − iK rot (0) | ε i (cid:105)| m (cid:105) , E = ε i + m Ω where H effrot | ε i (cid:105)| m (cid:105) = ε i | ε i (cid:105)| m (cid:105) and ˆ n | m (cid:105) = m | m (cid:105) H effrot commutes with ˆ n and (2) the operator e − i Ωˆ nt e − iK rot ( t ) e i Ωˆ nt is an analytic function of time t . While explicit expressions for K rot and H effrot can be obtained using the high frequencyexpansion, the above expressions should hold even when the HFE does not converge, so long as the stated assumptionsare satisfied.By Floquet’s theorem, the kick operator K rot ( t ) must be time periodic with frequency Ω. Hence it should berepresentable in a Fourier series in powers of e − i Ω mt . Because in the rotating frame, factors of e ± i Ω t are accompaniedby the corresponding qudit raising/lowering operator τ ± , we posit that terms in the Fourier series with exp( i Ω mt )should induce transitions between qudit states separated by (signed distance) m . The kick operator should thendecompose into K rot ( t ) = d − (cid:88) m =1 − m +( d − / (cid:88) n = − ( d − / e i Ω mt | n + m (cid:105)(cid:104) n | B n + mn + h.c. , where the operators B ij act only on the spins. In the HFE, one sees that B ij are imaginary and not necessarilyHermitian. We shall assume these properties still hold even when the HFE breaks down.We shall model the effect of the kick operator on only the qudit states by supposing that matrix elements of B ij between two delocalized spin states are random numbers, with possible dependence on i − j . In short, we propose thereplacement Tr S e − iK rot (0) ρe iK rot (0) −→ exp( − iK ) ρ qd exp( iK ) , where K is a d × d Hermitian random matrix whose upper triangular part (excluding the diagonal) looks like( K ) mn = ig exp (cid:32) − α (cid:18) m − nγ/ Ω (cid:19) (cid:33) R mn , for random R mn ∼ Normal( µ = 0 , σ = 1) and g, α >
0. The factor of γ/ Ω was inserted so that exp( − iK ) wouldbe the unit matrix in the decoupled and infinite frequency limits. The average over all realizations of K mimics thenonunitarity of the partial trace over the spins S .For example, we find good fits to the symmetrized occupations for the following values of the parameters, setting g = 1: L γ α
16 0 . / (cid:112) /
12 1 / . / . / (cid:112) /
12 1 / . / (cid:112) /
12 1 / ● ● ● ●■ ■ ■ ■ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ⨯ ⨯ ⨯ ⨯⨯ ⨯ ⨯ ⨯⨯ ⨯ ⨯ ⨯⨯ ⨯ ⨯ ⨯ n p n Symmetrized qudit occupationsScaled coupling γ = / L / ● L = (
10 realizations ) ■ L = (
10 realizations ) ◆ L = ( ) ▲ L = γ = ⨯ α - = ⨯ α - = ⨯ α - = ⨯ α - = We conjecture that the correct form in the limit of large L is( K ) mn = i exp (cid:32) − cL (cid:18) m − nγ/ Ω (cid:19) (cid:33) R mn , where cc