Long-time memory effects in a localizable central spin problem

We study the properties of the Nakajima-Zwanzig memory kernel for a qubit immersed in a many-body localized (i.e., disordered and interacting) bath. We argue that the memory kernel decays as a power law in both the localized and ergodic regimes, and show how this can be leveraged to extract $t\to\infty$ populations for the qubit from finite time ($J t \leq 10^2$) data in the thermalizing phase. This allows us to quantify how the long-time values of the populations approach the expected thermalized state as the bath approaches the thermodynamic limit. This approach should provide a good complement to state-of-the-art numerical methods, for which the long-time dynamics with large baths are impossible to simulate in this phase. Additionally, our numerics on finite baths reveal the possibility for unbounded exponential growth in the memory kernel, a phenomenon rooted in the appearance of exceptional points in the projected Liouvillian governing the reduced dynamics. In small systems amenable to exact numerics, we find that these pathologies may have some correlation with delocalization.


Introduction
Central spin models are ubiquitous in physical and chemical settings, from electrons with hyperfine coupling to nuclear spins inside quantum dots [1,2,3,4], to nitrogen-vacancy centers in diamond [5,6,7]. Depending on the couplings in these systems, the central spin may have long-lived, slow decaying dynamics suitable for quantum information applications. The role of the bath in these cases is relegated to modeling decoherence, and has not traditionally been considered to be important. The bath is usually taken to be non-interacting, an assumption which has proven fruitful in the development of analytical [8,9,10,11,12,13,14] and numerical [15,16,17,18,19,20] techniques. As such, these classes of baths-whether composed of bosons [21] or spins [22]-are by now reasonably well understood [22,23,24,25,26,27,28,29,30].
Recent research has brought new focus to modifications of the bath by adding, for example, intra-bath interactions and disorder. With these additions, the bath alone can exhibit novel dynamical phases such as many-body localization (MBL), which serves as a basis for nonergodicity in generic systems with strong disorder. Upon coupling to a bath, the long-ranged mediated interactions between constituents of the bath can push the bath towards delocalization. Recent work [31,32] has shown that a single qubit coupled centrally to a 1D MBL spin chain can preserve localization provided that the magnitude of the central coupling decays fast enough with the size of the bath. Delocalization can be achieved by a sufficiently strong magnitude of central coupling, which the authors of [31] took to be signaled by quantum chaotic energy level statistics. However, it was noted in [32] that the nature of the delocalized phase is unclear, as it could be nonergodic. This could be reflected in the long-time value of the central qubit's population not reaching the thermal expectation but, as was found in [10] studying integrable central spin models perturbed away from integrability, limitations of bath size prevent a definitive conclusion. An impediment is that at strong couplings, analytical and numerical approaches become scant due to the presence of interactions in the bath and to the star-like geometry of problem.
The added complexity has a drawback in that such systems quickly become intractable computationally, even for small bath sizes of ∼ O(30) degrees of freedom. This is due to the exponential increase of states in the Hilbert space that are involved in the dynamics. Moreover, large intra-bath interactions and disorder can radically change the timescales of the bath and invalidate perturbative approaches to bath dynamics. In this context, the case of MBL (in one dimension) is special in that it allows for a non-perturbative description in terms of 'l-bits' [33,34,35,36]. Owing to this and to the slow growth of entanglement entropy [37,38], dynamics in the localized phase of MBL systems are by now well explored numerically and analytically [39,40,41,36,42]; however, these approaches generally fail on the ergodic side of the transition.
The dynamics of extended, thermalizing many-body systems are typically very difficult to simulate exactly due to the rapid growth of entanglement. This is, for example, the limiting factor in methods based on a tensor network ansatz for the wavefunction in which the bond dimension bounds the amount of entanglement entropy that can be captured. A reasonable strategy then would be to extend the timescale of the converged simulation using information that can be computed on the timescales before the breakdown of the numerical method. Such an approach had been used successfully in the past to find the steady state behavior of quantum impurity systems [43,44,45], and to show the existence of bistability in the Anderson-Holstein model [46]. In those applications, the nontrivial dynamics of the impurity could be described exactly using a memory kernel, derived using the projection operator formalism described by Nakajima, Zwanzig and Mori [47,48,49].
While the Nakajima-Zwanzig theory is formally exact, it is oftentimes more demanding than other formalisms to describe the dynamics because of the time-nonlocal memory kernel that naturally arises in their approach. It only becomes computationally useful if the nonlocality can be restricted, e.g. large timescale separation between bath and system dynamics lending to Markovian approximations, or if memory kernel decays sufficiently rapidly such that it can be truncated for times ≥ t c , where t c is the cutoff time.
The use of memory kernels to study dynamics in central spin systems has seen various degrees of success [24,50,26,13]. For analytical tractability, such studies are usually restricted to noninteracting baths without disorder and the memory kernel is expanded perturbatively. In the cases where such expansions are valid, it has been found that the memory exhibits nonexponential decay at long times, with long-time averaged population consistent with a nonergodic dynamics [24]. However within the perturbative approach it is found that at higher orders of the expansion, the memory kernel can display secular (unbounded) growth [24,26]. In this work, we shall go beyond these approaches, taking into account the presence of bath-bath interactions along with random disorder and directly computing the memory kernel, therefore bypassing the possibility of pathological behaviors in the perturbation.
In this paper we study the memory kernel of a two-level system immersed in a bath modeled by a many-body localizable spin chain. We do so with two goals in mind: to assess the feasibility of extending the system dynamics from short time calculations when analytical and direct numerical approaches to compute the system dynamics fail (i.e. on the thermalizing side of the MBL transition); and to understand how interactions and disorder in the bath affect the memory kernel in properties such as timescales and tail behavior.
To this end, we will work with a previously studied model [31] of a qubit (τ x,y,z ) coupled to a disordered Heisenberg chain of L spins-1/2 (σ x,y,z i ): where theτ andσ are Pauli matrices. The bath HamiltonianĤ B corresponds to the disordered, isotropic Heisenberg chain, where we take J = 1. The system-bath coupling termsV are likewise given by the Heisenberg interaction, with magnitude scaling as γ/L to ensure that localization can occur for finite γ. We shall refer to γ as the strength of the central coupling. The random longitudinal fields h i are drawn independently and uniformly from [−W, W ]. The data we present here will be restricted to W/J = 6, chosen such that the bath is localized for γ = 0 and experiences a central couplinginduced delocalization [31] around γ ≈ 5. Finally, the magnetic field is set to Ω = 0.
Thus the qubit has no intrinsic dynamics and is instead entirely dependent on the magnitude of the Overhauser field it experiences from the bath. As noted in [31], the interacting central spin problem of Eq. (1) can be realized in dipolar spin ensembles. Other platforms that allow for experimental realizations of this model include programmable quantum simulators [51], NMR experiments with triphenylphosphine [52], or in superconducting qubit circuits [53]; these approaches offer a high degree of control, for example by enabling the control of intrabath interactions, random disordered Zeeman fields, and the strength of coupling γ.
The structure of this paper is as follows: we shall first define the memory kernel for reduced dynamics and consider the role of disorder averaging; then we shall analyze the physics underlying the memory kernel at short, intermediate, and long times; and finally we shall discuss the potential for the memory to be used to augment short-time experimental or numerical data.

The Nakajima-Zwanzig equation
We quickly review the basics of the projection operator approach to generalized quantum master equations. Any given Hamiltonian can be split into contributionsĤ S acting only on the system,Ĥ B acting only on the bath, andV coupling the two. We will use the term "bath" as a shorthand for the set of physical degrees of freedom surrounding the central qubit. In particular, we do not assert the character of the bath to be unchanged by coupling to the system. To each of the three aforementioned operators is associated a corresponding Liouvillian superoperator ( generating dynamics for the density matrix Oftentimes one is interested only in the dynamics of the system, in which case the bath degrees of freedom can be projected out by tracing over the bath on both sides of the equation, where the bath trace is This is used to define the system reduced density matrix, We shall additionally assume that the initial state is factorized, i.e.ρ 0 =ρ S,0 ⊗ρ B . By taking the bath trace defined in (3) on both sides of (2) and using Tr B L B = 0, we arrive at the exact expression which is an equation of motion forρ S (t) that explicitly depends on knowledge of the time evolution of the full system and bath. This equation of motion can be closed, i.e. involving onlyρ S (t), by using Dyson's identity (see [54,55]): The memory kernel superoperator is formally defined as In the above equation, the projection superoperator is taken to be P· ≡ Tr B ( · ) ⊗ ρ B and Q = I − P is its complement. It is useful to define the system reduced propagator (superoperator) such that Knowledge of U S allows for the generation ofρ S (t), and lets us write a Nakajima-Zwanzig equation [56] involving only objects of one type, i.e. superoperators: In this form, it becomes clear that one can solve for K directly from U S . Note that no approximations have been made and the dynamics generated by solving (9) and (7) are equivalent to solving (2) with the stated assumptions on initial conditions. The derivation however, benefits from a simplification made possible by the form of the model Hamiltonian in (1). Bath traces over the interaction Liouvillian L V with respect to a bath stateρ B of fixed magnetization will be zero due to the conservation of total magnetization in the model, and if we chooseρ B to have zero magnetization. Therefore the validity of (7) is not restricted to solely thermal baths (ρ B ∝ e −βĤ B ) nor bath eigenstates ( ρ B ,Ĥ B = 0).
The memory kernel K and the system propagator U S , being linear mappings from the system Hilbert space H S to itself, can be represented as (dimH S ) 2 × (dimH S ) 2 matrices. Requirements on unitarity and hermiticity, along with the decoupling of populations and coherences in this magnetization-conserving model, means that U S is described by only two independent entries when the focus is solely on population dynamics. The same extends to K by virtue of its relation to U S . The two entries of U S are computed by two independent instances of the initial system stateρ S (0): one from the population of the |0 state whenρ S (0) = |0 0|, and the other from the population of the |1 state whenρ S (0) = |1 1|. The initial bath state is the same in both cases, with definite magnetization M B = 0. Because the total magnetizationM z =τ z + iσ z i is conserved, these two trajectories must reside in independent parts of Hilbert space. They are then combined in solving for the memory kernel, which can be done in the time domain by discretizing the integro-differential equation (see the supplementary materials for details). This, while posing no problem for the projection operator formalism, leads to a strange scenario where the central qubit dynamics restricted to one symmetry sector will depend on information from another, disjoint symmetry sector.
To skirt around this unsavory philosophical scenario, we can focus on only the population of the |0 state of the central qubit. Using the projection operator Pρ = (|0 0| ⊗ρ B ) Tr[(|0 0| ⊗Î B )ρ], one can repeat the same steps as before and obtain the scalar memory kernel for a single disorder realization as satisfying the integro-differential equation or its Laplace-domain equivalent Focusing on the population p 0 (t) of single state allows us to work with a scalar memory kernel K(t) and simplifies the calculations. We will focus exclusively on the scalar memory kernel for the remainder of this paper. While this may be an unconventional choice of projector and therefore also of a memory kernel, we stress that the Nakajima-Zwanzig equation in its most general form does not depend on the choice of the P. The only requirement is that the same observables of interest are contained in the domains of the different projectors. ‡ Note that the memory kernel is akin to the self-energy for the reduced density matrix. Solving for it is then tantamount to solving the exact problem. Yet there are still advantages to working with the memory. For one, because of its relationship with the central qubit's populations it is in principle a measurable quantity. There is also the possibility for the memory to decay on timescales different from that of the populations. Should the memory decay much faster, then it may be possible to leverage the timescale separation to reduce the computational effort required to solve for the system dynamics at longer times.

Disorder averaged memory
Given that we are interested in disordered systems, suitable definitions of a memory kernel associated with different disorder realizations depends on the quantity of ‡ See [57] for a detailed demonstration of the equivalence of dynamics generated by different forms of generalized master equations resulting from the interplay of projections and the presence of conserved quantities. experimental interest. The difference depends on when the disorder averaging is performed. We denote by K avg the case where the population p 0 is averaged over the disorder (p 0 ) before solving for the memory kernel, satisfying The other case, where the memory for disorder realization is found and then averaged, is denoted by K. This latter case is relevant should one decide that the observable of interest is the memory kernel itself, which is in principle possible since it is directly computable from the populations. It is not a priori clear how these two definitions are related. A reasonable guess might be that, upon disorder averaging, the two definitions are equivalent. We argue that this is not necessarily correct. Suppose that for every L the disorder-averaged population p 0 (t) exists, with initial condition p 0 (0) = 1. The trajectory of the population for a single instance of disorder will have deviations from this average value, p 0 (t) = p 0 (t) + δp(t). Since populations must be positive at all times, so should their Laplace transforms for real, positive z. Using (12), the positivity of the Laplace transforms allows us to write Since the exponential function is entire, the term in parentheses can be expanded as a series, Averaging this expression over disorder, we will have the n = 1 term vanish by definition of δp. But all higher order terms-particularly ones with even powers-are not guaranteed to vanish. The consequence is that K avg = K for finite L.
The situation is modified in the thermodynamic limit owing to self-averaging. Intuitively, a small subsystem interacting randomly with N ≫ 1 degrees of freedom should have deviations from its mean behavior that decrease as N increases. As a result, when the environment is sufficiently large, a single realization of the random interaction should typically yield results close to the mean. This statement was recently demonstrated [58], showing that the system reduced density matrix enjoys the typicality property for system-bath interactions modeled by certain classes of random matrices. Importantly, [58] showed that this self-averaging property holds at least up to a timescale T that increases with L. Thus if p 0 (0 ≤ t ≤ T ) is self-averaging, so must K(t) on the same interval, since to solve for K(t) up to time T in (11) requires only p 0 (t) on [0, T ]. In figure 1b we show the root-mean-squared fluctuations of p 0 (t) deeply in the thermalizing phase of the bath-disordered Hamiltonian (1), and observe that they indeed decrease with increasing bath size. Extrapolating to the thermodynamic limit, we should therefore have self-averaging of the reduced density matrix of the central qubit. Then by extension the memory must self-average too. This can be seen from (15), where fluctuations of a single realization of K(t) has deviations from K avg (t) that are bounded by the magnitude of the fluctuations in the population δp(t) = p 0 (t) − p 0 (t). In figure 1a, we find that K(t) and K avg (t) generally tend to differ by |K(t) − K avg (t)| ∼ O(10 −3 (γ 2 /4L)) up to timescales t O(10 2 ) for the system sizes we can simulate. We observe that this deviation can diverge exponentially with a finite number of disorder realizations at long enough times, a phenomenon which we will return to in section 2.3. Barring that, the self-averageness of the population p 0 (t)which yields K avg (t) = K(t) =⇒ K avg (t) = K(t) in the thermodynamic limit-gives us an alternate window into understanding how the memory kernel behaves. For the remainder of this paper, we shall mostly discuss K avg (t) as we are interested also in the dynamics of the averaged population.

Results
We implement time evolution by approximating e −iĤt with Chebyshev polynomials [59,60]. To reduce computational costs, we use the conservation of total magnetization M z =τ z + iσ z i in the model, allowing us to restrict the dynamics to the symmetry sector withM z = −1. The system is prepared in theρ S,0 = |0 0| state, while the bath stateρ B is initialized to be a Neel state, | · · · ↓↑↓↑↓ · · · . We expect similar results should we choose different initial states within the sector ofM z = −1.
A (matrix) memory kernel K with n independent entries can be computed directly from the populations [56] using n different initial conditions, for each disorder realization. In this sense, there is added computational benefit to restricting our discussion to only the scalar memory kernel K(t).

Short times
We can leverage the self-averaging property to gain some understanding of the short time behavior (figure 1c) of K avg (t). The derivatives of K(t) at t = 0 for a single disorder realization can be found straightforwardly (see the supplementary materials) from those of p 0 (t), with the lowest orders being where f (n) denotes the n-th derivative. After averaging over disorder with an initial Neel state in the bath, we have where J = 1 in our model, and L is the number of spins in the bath. We see that the disorder strength W sets the initial decay rate 1/τ K . This can be roughly estimated for large W from Fermi's Golden Rule, using our argument that disorderaveraging effectively gives us a continuous spectrum with an effective (root-meansquared) bandwith ∼ O(W √ L), and a coupling strength ∼ (γ/2L) 2 . While W sets the decay timescale for K avg (t), it is the quantity γ 2 /4L that sets the overall magnitude of K avg (t) and so dictates the timescale for p 0 (t). We expect so from the following scaling argument: Assume that the memory kernel converges to a limiting form in the thermodynamic limit as where k(t) is independent of L and has a short time expansion given by (17). From the Nakajima-Zwanzig equation, we rescale the time to t ′ = (γ r /L s )t and obtain dp ′ where p ′ 0 (t ′ ) ≡ p 0 (t ′ /(γ r L −s )). We seek exponents r > 0 and s > 0 such that p ′ 0 (t ′ ) will vary on the timescale ∆t ′ ∼ 1. With the rescaled time, the k(L s τ ′ /γ r ) appearing in (20) will have largely decayed by τ ′ ∼ γ r L −s /(W/ √ 3), a timescale much faster than that of p ′ 0 (t ′ ). Hence we can approximate p ′ 0 (t ′ − τ ′ ) in (20) as a constant, and estimate the strength of memory effects by integrating k(L s τ ′ /γ r ) up to its decay time. This is roughly given by We require s = 1 in order to have a converged p ′ 0 on the timescale of t ′ in the thermodynamic limit. Furthermore, r = 1 so that a trivial rescaling of the Hamiltonian H → αĤ would not alter the strength of the memory term. Thus we argue that the dynamics of the central qubit should proceed on the timescale τ p 0 ∼ L/γ, consistent with our initial assumption that the population dynamics proceed much more slowly than does its associated memory kernel. We show this rescaling of time in figure 1b and in the inset of figure 1d, where the former shows the fluctuations of p 0 (t) between disorder realizations for different system sizes at fixed γ = 10, and the latter shows p 0 (t) for fixed L = 12 across γ. These figures show that the lowest moments of the populations align on the timescale τ p 0 ∼ L/γ. This result is also consistent with the result of [31] on the central qubit's autocorrelation function, dτ τ z (t + τ ) τ z (τ ) , where it was observed that there is an accumulation of spectral weight near ω ∼ γ/L.
With a clear separation between τ K and τ p 0 , one may wonder whether the central qubit can be described by an effective master equation. At least deep in the localized phase, the bath is too slow to act as an effective reservoir for the central system. Correlation functions of the bath are argued [61,62] to decay as a power law t −ζ with 0 < ζ < 1, which makes memory effects crucial in dictating the behavior of p 0 (t) at long times. We will return to discuss the long time behavior of the memory kernel below in section 2.3.

Intermediate times
As seen in figure 1d, the memory past Jt 1 takes on different behaviors depending on the coupling strength, with increasingly damped oscillations as the combined system and bath transitions from localization to thermalization. The inset of figure 1d shows that this behavior is not observable when looking solely at the populations. The oscillation is dominated by frequencies in the range ω ∈ (4, 6), close to the disorder strength W = 6. Such oscillations are not a feature unique to an interacting bath. They show up in the non-interacting limit J = 0, in which the memory to lowest order in γ can be approximated by where γ ⊥ = γ/2. We see that oscillations are linked to the finite bandwidth W of frequencies in the bath [43], which arises from precession about the local field on each site, (h i /2)σ z i , and h i ∈ [−W, W ]. When interactions in the bath are turned on, we would expect them to provide a small renormalization to the precession frequencies, as we are working with a hierarchy of scales such that W ≫ J > γ. This assumes, of course, that the bath dynamics are approximately describable with a precession picture even in the presence of bath interactions.
To justify this picture more formally, we can leverage the description of MBL systems in terms of quasi-local integrals of motion, which form the effective bath degrees of freedom that exhibit precession. At intermediate times and at weak coupling, the memory kernel can be approximated by bath correlation functions [63,64], In the MBL phase, the bath spin operatorsσ ± i have large overlaps [33] with quasi-local operatorsΘ x,y,z i with which the bath Hamiltonian can be written as [33,35,62] where the operatorsΘ x,y,z i follow the Pauli commutation relations. The bath correlation functions oscillate according to ε i , at least when the bath is strongly localized. The distribution of ε i will therefore dictate the intermediate-time behavior of the memory kernel. For instance, if the distribution has sharp cutoffs like in the case of box disorder, then it can be expected that the memory will display oscillatory behavior whenever the stated approximations are applicable. We note that the picture of precessions is complicated at later times by dephasing mechanisms arising from interactions-the 2body J i,j terms and higher-in the bath. Therefore, measurement of the memory kernel will yield some information on the parameters entering the bath Hamiltonian (24).
At the other extreme, where the system strongly couples (γ 5 for the value of W = 6 we have shown in figure 1d) to the bath, the localization assumed above breaks down [32]. That is, the bath interactions mediated by the qubit are strong enough that the bath cannot remain "close" to its initial state, so the expansion of the memory in terms of bath correlation functions no longer holds. In all, the contribution of the bath to the system dynamics can no longer be parsed into contributions from (nearly) independent oscillators. Instead, delocalization evidently serves to homogenize the influence of the bath, smoothing over the randomness from the local fields h i , and damping out oscillations in K avg (t) as observed in the red curves of figure 1d.

Long times
Within the particular parameters we have chosen to study in this model, we define "long times" to correspond to Jt 10, a time past which the coherent oscillations in the bath have dephased. For the purpose of extrapolating the dynamics, it is crucial to understand how quickly K avg (t) decays, if it even does so at all. However, since we can only numerically average over a finite number N of disorder realizations, we in the tail portion K avg (50 ≤ t ≤ 100) decays as 1/ √ N , and moreover decreases with increasing system size as would be expected from self-averaging systems. In the above equation, we use the notation g(t) [ The persistence of the finite N noise makes it difficult to conclusively show numerically whether K(t) decays as algebraically or exponentially. While in section 2.1 we argued for a power law decay for the weakly coupled, localized phase based on known phenomenology of MBL, this approach cannot work for the strongly coupled, thermalizing phase. In the absence of weak coupling perturbative expansions we now turn to the self-averaging relations K avg ∼ K ∼ K to attempt to extract insights about the thermalizing phase. Doing so requires discussion about the memory kernel for a single realization of disorder, which is what we shall focus on for the remainder of this subsection.
For certain realizations of {h i }, we observe an increasing likelihood for the memoryboth scalar-and matrix-valued versions-to display unbounded exponential divergences with increasing coupling γ. We can verify the divergence for small system sizes L 6, where the Laplace transformed memory kernel can be computed directly to yield the memory as a sum over simple poles, some of which with positive real parts. Such contributions-which are necessary in order to correctly reproduce the population dynamics-lead to an unbounded exponential increase of the memory for particular values of the coupling and magnitude of disordered fields. We will return to discuss the origins and implications of such pathological behavior in section 4.
We can motivate the consequences of exponentially growing contributions to K(t) by examining the structure of the poles of its Laplace transform, K(z). Because the Hamiltonian is real and Hermitian, poles of K(z) are given by a real polynomial (see the supplementary materials for details). The polynomial will only involve terms of even powers, z 2n , because p 0 (t) = p 0 (−t). Thus if a pole s n exists with residue r n such that Re s n = 0, it must be the case that poles −s n , s * n and −s * n must exist with residues r n , r * n , and r * n respectively. Based on the distribution of the pole structure, any exponentially dampened part of the memory (Re s n < 0) must be accompanied by an exponentially growing counterpart. We posit that in the thermodynamic limit one of two situations must hold: 1) all off-axis poles converge towards the Im z axis as L → ∞, or 2) some poles still exist off-axis, which because of the conjugate pairs, contributes both exponential decay and growth. In the first scenario, there are no isolated poles to cause exponential decay. In the second scenario, any exponential decay is masked by exponential growth. Moreover, even if Re s n > 0 poles cancel upon disorder averaging, the same would happen to the Re s n < 0 poles by virtue of the relationship between residues discussed above. Therefore we argue that even in the thermalizing phase, the memory kernel for the dynamics we have defined should not exhibit exponential decay in the limit as L → ∞. This leaves open the possibility of power-law or stretchedexponential behavior. In the next section, we will use infinite-time data from exact diagonalization to show that the long-time behavior of the memory is consistent with a power-law decay. Finally, we reiterate the importance of the order of limits in this problem. They must be taken as lim t→∞ lim L→∞ lim N →∞ (26) to ensure that, reading from right to left, the population-and therefore the memory kernel-does not recur and to ensure the validity of the approximation K ≈ K avg .

Extracting long time information from the memory kernel
The memory kernel has a direct relation to steady state values of the reduced density matrix, provided that a steady state exists [43,44,46,45]. While the past work was done using all d 2 ×d 2 elements of the (matrix) memory kernel, we can import their ideas to the scalar memory kernel and a single element of the reduced density matrix. From the relationship between p 0 (z) and K avg (z), we can use the final value theorem to find lim z→0 z p 0 (z) = lim If κ(t) decays sufficiently quickly, we can extrapolate the z → 0 limit from the finite times accessible from numerics. However, as we argued in the previous section, the A long time tail of K avg (t) would have non-negligible contributions to the dynamics, and therefore much care has to be taken in its use for extrapolations. In lieu of a cutoff approximation, we turn again to the definition of K avg , We take an ansatz for the memory at small z, where 0 < ζ < 1 and the long time limit of the average population p 0 (t) shall be denoted as p ∞ . Note that we had argued in the previous section at least for the absence of exponential decay of the memory kernel in the thermodynamic limit, based on the structure of the poles in Laplace space. The presence of terms like z ζ is consistent with long-time behavior as To extrapolate the long time populations, we compute κ(t) defined in (27) and approximate its Laplace transform This result is then fitted using (30) to find p ∞ and b 0 and the amplitudes a n . Such an approximation for the Laplace transform is admissible only if κ(t) has decayed to sufficiently small values at t = t max , and for z t −1 max . We find that the results of using such an extrapolation procedure agree well with the values from independent calculations using exact diagonalization (figure 3a). Thus we are able to obtain estimates for the long-time population of the central qubit for system sizes (L 16) larger than those obtainable through exact diagonalization. In particular, this allows us to see how the central qubit approaches the thermalized limit p 0 = 1/2 with increasing bath size. In figure 3c, p ∞ is consistent with power law decay p ∞ − 1/2 ∼ L −1.03 , which is in line with the scaling given by the infinite temperature phase space average, measuring the relative sizes of the Hilbert spaces for eigenstates occupying |0 and |1 . We stress that because the memory must decay with time, this procedure cannot be used in finite systems for a single disorder realization, as the population will generally not reach a steady state in such circumstances. Furthermore, since we have estimates of the true value of p ∞ obtained independently from exact diagonalization, we can compare (30) to a more generic alternative where κ(z) is analytic about z = 0. Such is the case if κ(t) were to, for example, decay exponentially or faster. For different system sizes and disorder distributions, we have found that only the power-law ansatz is able to smoothly interpolate between known z = 0 values of κ(z) from exact diagonalization and z > 0 values of κ(z) calculated from finite time dynamics (see the supplementary materials for an example). Thus, while we have been unable to mathematically prove the existence of a long-tail in K avg (t), we have at least found numerical corroboration for the validity of our claim.
We note that, at least for L ≤ 14, we find that the exponent ζ is system size dependent, for both box (figure 3b) and Gaussian distributed disorder. In the absence of intrabath interactions (J = 0), we observe that the integrated kernel κ(t) acquires a large oscillatory component with a decaying envelope at long times for γ = 10, which dominates over the κ(t) ∼ t −1−ζ behavior seen with J = 1 (cf. figure 3d), see section four of the supplementary materials. We further argue in the supplementary materials that if one takes the bath to initially be at infinite temperature, there will be a temporal powerlaw decay of the memory as ∼ t −3 which implies that ζ → 1 in this limit. Altogether, this suggests that ζ is at least a quantity dependent on intrabath interactions as well as the initial state; we cannot clarify whether there is a limiting value as L → ∞ for initial states of fixed energy density, such as that considered in this work.
One may wonder what advantage this method confers to obtaining infinite-time populations, compared to simply simulating the population dynamics to longer time. For one, it is not always clear the timescales at which one can be sure that the system will have relaxed. This point is made more salient by the possibility of small, longtailed memories which implies similar behaviors in the population dynamics. In this  [31]. The dashed black line is the asymptotic behavior of the boundary as argued in Ref. [32].
work, we have argued that it suffices to be able to observe whether the memory has reached the regime of power-law decay, at which point one can use (30). We stress that power-law behavior may become more apparent at earlier times in the memory than compared to the population, such as what we have observed in this work. While the t −ζ−1 contribution to κ(t) may be subtle-on the order of 10 −3 in all system sizes and disorder distributions we examined (see figure 3d for an example)-it is always possible to systematically improve its resolution simply by performing more disorder averaging.

Unbounded exponential growth of the memory kernel
We return now to the observation made in section 2.3 about memory kernels growing exponentially in time for certain realizations of the disorder. As seen in figure 4a, this can show up in the disorder averaged memory K(t), which can only be approximated via sampling over a finite number of disorder realizations. In figure 4b we show the maximum real part of the poles-corresponding to the maximum rate of exponential growth ν-for specific set of {h i } with L = 4. Intriguingly, ν is not monotonic with respect to γ, and displays square root singularities when going from ν = 0 to finite ν.
The sharpness of these singularities even with L = 4 indicates that they should not be associated with thermodynamic phase transitions. Instead, we believe they stem from exceptional points (EPs) in the generator of projected dynamics, QLQ, which are related to generalized avoided crossings. This generator is responsible for the time evolution of the memory kernel, as seen in (7). By choosing to focus on only a subset of all the physical degrees of freedom in the problem, we were forced to define projection operators P that are not self-adjoint in the space of operators [65,66]. For example, in operator space the projection operator associated with the scalar memory kernel is P = ||0 0| ⊗ρ B )(|0 0| ⊗Î B |, where the adjoint of the operator state vector has action (Â|B) = Tr Â †B .
The condition of being self-adjoint Liouville space is Writing P in this way, it is clear that even if we project on to a thermal state of the bath,ρ B ∝ e −βH B , the projector P still cannot be self-adjoint unless the bath is in an infinite temperature state. Thus the projected Liouvillian QLQ is also not self-adjoint, a property which allows EPs to occur. We have verified that the same phenomenon occurs even if we work with larger projection superoperators leading to matrix-valued memory kernels. We have additionally verified numerically that features unique to EPs such as the coalescence of eigenvalues and self-orthogonality are also present (see supplementary materials). Interestingly, we note that the region in (W, γ)-space (figure 4c) for which MBL is predicted to be stable in the thermodynamic limit appears to be correlated with a suppressed ν. While we are currently unable to prove that this is not a coincidencesuch system sizes cannot inform us about the stability of MBL -it is possible that this provides a window into the character of the eigenstates, which are argued to be radically altered at large enough γ due to percolating networks of resonance states [32]. At the same time, it is known that the presence of exceptional points limits the radius of convergence for perturbative expansions [67,68], and is postulated to be linked to quantum phase transitions [69,70]. To fully explore any link between exceptional points, delocalization, and the breakdown of perturbative approaches to MBL will require a separate, in-depth study.
Heuristically speaking, delocalization with increasing coupling is the result of singular behavior in the full Hamiltonian, a fact which should be reflected in both the eigenstates and the spectrum. In finite systems, these may be isolated occurrences whose singular properties are smoothed out upon taking expectation values. Our numerical observations suggest that the non-Hermiticity of the projected Liouvillian is highly sensitive to such singularities. We suspect this may be further indication of a deeper connection between localization and long-time pathologies in the memory kernel, but we are unable to clarify the underlying physics at this time. However, we will note that the situation may be altered by introducing a large bias on the central qubit, e.g. Ωτ z , the analysis of which we will leave for future work.

Discussion and conclusions
In this work we have undertaken the study of the time-nonlocal memory kernel describing how a many-body localizable "bath" affects the population dynamics of a central qubit. While the memory is formally defined in terms of Liouvillians, the dimensions of which quickly grow to be computationally intractable with increasing size of the Hilbert space, we are able to compute it numerically exactly from existing methods for simulating dynamics in closed quantum systems [56]. We note in passing that the method we use in this work is general, and can easily be formulated to describe the dynamics of the central qubit's coherence, as might be relevant for some recent NMR experiments [52]. With this method we are able to directly examine the behavior of the memory kernel, parsing it into three regimes: short, intermediate, and long times.
On short timescales (Jt 1) is where the majority of the memory's decay occurs, irrespective of whether localization (at small γ) or delocalization (at large γ) is present. Properties of the memory on this timescale largely dictate the timescale of the dynamics for the central qubit's populations. We note that (17) holds for arbitrary unbiased (i.e. zero mean) distributions of onsite fields with variance W 2 /3 and bath-bath interaction strength J. On intermediate timescales (Jt 10) in the localized phase, the memory should exhibit dynamical signatures that result from the distribution of effective couplings for the emergent local integrals of motion describing the localized bath. For example, if the disorder distribution has sharp cutoffs, then this is manifest as oscillations in the memory. These oscillations are damped out as γ is increased, tuning the system and bath into the thermalized phase. This behavior strongly depends on the distribution of disorder, as well as on the presence of bath-bath interactions. Finally, at long times (Jt 10) we observe pathological exponential divergence of the memory kernel for certain realizations of disorder, deep in the thermalizing phase. We find that this comes from exceptional points in the projected Liouvillian generating the dynamics of the memory kernel, which come about at real values of the coupling γ due to the non-Hermiticity of the projection superoperator used to define the projected dynamics in the Nakajima-Zwanzig formalism. Unlike past work [65] that treated such exponential divergences as unphysical and should therefore be discarded, we have taken the view here that the divergences have a meaningful impact on the population dynamics. We argued that after disorder averaging the memory kernel, such pathological behaviors should preclude any exponential decay of the memory. Instead, we find that the tail of the memory is consistent with a power-law decay ∼ t −2−ζ , where 0 < ζ < 1. We find that this form still holds true for different distributions of disorder. However, in the noninteracting bath case of J = 0, the strictly power law decay appears to be replaced with an oscillatory component with a decaying amplitude that we find to be consistent with a power-law. In the interacting (J = 1) case, such a power-law ansatz allows us to extract estimates of the disorder-averaged infinite-time population of the central qubit, solely from finite-time simulations. While such a procedure was shown in the past to work well when one could define a cutoff time for the memory kernel [43], here we have argued for the possibility that no cutoff time exists and demonstrated a proof-of-concept approach for extracting the infinite time populations in such a scenario.
In the model we have studied in this paper, we have taken the central coupling to scale to zero as γ/L, in accordance with Refs. [32,31] which have argued for its necessity to perturbatively preserve localized eigenstates. As a consequence, we have argued that there arises a separation of timescales between the population dynamics (τ p 0 ) and its associated memory kernel (τ K ). Should we repeat our arguments from section 2.1 with a central coupling scaling as γ/L q , we find that these two timescales remain separated for q > 1/2, but coincide for 0 < q ≤ 1/2. It is not clear whether such a separation of timescales-where τ p 0 ≫ τ K as L → ∞-is required for the preservation of localization. Heuristically speaking however, having τ p 0 ≫ τ K does not appear at first glance to be strong enough to preserve all aspects of MBL. One of the dynamical hallmarks of MBL is a logarithmically slow spreading of entanglement, i.e. spins on sites i and i+L/2 become entangled after a timescale ∼ exp(L/2ξ) with ξ being the localization length [36]. Based on our view of the system dynamics from the memory kernel, the interaction between these two sites mediated by the central qubit should proceed on a timescales growing as a power of L, which is much shorter than the dephasing time ∼ exp(L/2ξ) and thus may accelerate the dephasing process responsible for the slow dynamics in the MBL phase. However, it was noted in Ref. [31] that the central qubit at best facilitates a subextensive transport of magnetization which augments, but does not destroy, the logarithmic growth of bipartite entanglement.
Our work also raises tantalizing questions about possible connections between poles of the Laplace-transformed memory kernel and thermalization/delocalization. To this end, some work [67,71,72,73] has been done to connect the proliferation of exceptional points in non-Hermitian systems to the appearance of quantum phase transitions and chaos. By focusing on a subpart of a closed system, we are forced to consider non-Hermitian Liouvillians giving rise to exceptional points in the space of operators. Explorations in this direction may benefit from insights from the physics of Feshbach resonances. Of course, we are severely limited by the system sizes amenable to numerical studies, thus we are able to do little more than remark on the coincidences we observe.
On the more practical side, we have demonstrated that there may be enough information from finite time dynamics to yield knowledge about long time limits, should they exist. While we have only demonstrated the extrapolation to t = ∞ of the population of the central qubit, we should in principle be able to use the same memory kernel and the Nakajima-Zwanzig equation in (10) to extend the computed dynamics to longer times. That this is even possible should not be too surprising, given that (10) when discretized over time gives the same form as the ansatz underlying linear prediction [74,75], a method widely used for extending dynamical calculations. What we have shown in this work is that there may be more physical content in such a procedure than was previously appreciated. To explore these ideas more thoroughly warrants careful attention, particularly in regard to stability and applicability, which we shall leave for future work.
Finally, we note that any possibility of a pathological memory kernel at real γ can be erased by choosing to work with self-adjoint projection superoperators P. One may be interested in doing so, for example, in order to approximate system dynamics from low order, analytical expansions of the memory kernel. In that case it would be beneficial to know that the error introduced by the approximation is not exponentially divergent with time. It is as yet unclear whether self-adjoint projectors necessarily yield improvements, since pathological behaviors can still occur for complex couplings γ to limit convergence of naïve series expansions. We note, however, that previous work [76,26] saw benefits from applying symmetry-adapted "correlated projectors"-which, we should point out, are manifestly self-adjoint in Liouville space-to low order expansions of the memory kernel. We leave clarification of this point for future work.

Acknowledgments
We are grateful to Amikam Levy, Sebastian Wenderoth, and Michael Thoss for useful discussions. This research used resources of the National Energy Research Scientific Computing Center, a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.