Compelling Bounds on Equilibration Times -- the Issue with Fermi's Golden Rule

Putting a general, physically relevant upper bound on equilibration times in closed quantum systems is a recently much pursued endeavor. In PRX, 7, 031027 (2017) Garc\'{\i}a-Pintos et al. suggest such a bound. We point out that the general assumptions which allow for an actual estimation of this bound are violated in cases in which Fermi's Golden Rule and related open quantum system theories apply. To probe the range of applicability of Fermi's Golden Rule for systems of the type addressed in the above work, we numerically solve the corresponding Schr\"odinger equation for some finite spin systems comprising up to 25 spins. These calculations shed light on the breakdown of standard quantum master equations in the"superweak"coupling limit, which occurs for finite sized baths.


I. INTRODUCTION
The last decades have seen a major progess in the field of equilibration in closed quantum systems [1]. Concepts like typicality [2][3][4] and the eigenstate thermalization hypothesis [5,6] have been brought forth. Furthermore, it has been established that for an initial state ρ populating many energy levels, expecation values A(t) will generically be very close to their temporal averages for most times within the interval to which the average refers ("equilibration on average"). While the conditions for this statement to be true are rather mild concerning the observable A and the Hamiltonian H [7][8][9], the respective time interval may be very large. For specific observables it may, e.g., scale with the dimension of the relevant Hilbert space [10,11]. Moreover, concrete examples are known in which the corresponding equilibration times for physically relevant observables scale as T eq ∝ N α , α ≥ 1/2, where N is the size of the system. This result has been found for systems featuring long range [12] as well as short range interactions [13], albeit for a somewhat different definitions of equilibration times. In fact, already for mesoscopic many-body systems with standard interaction strengths, the required equilibration interval may be on the order of the age of the universe [12]. Thus, although the above statements in some sense establish equilibration under moderate conditions in the very long run, it is unclear whether or not this equilibration will ever occur in a physically relevant period of time. Hence, the question of an upper bound on this relaxation timescale has recently been much discussed. Since it is always possible to find mathematically well defined, permissible initial states that fully exhaust the above, unsatisfactorily large time interval, most contributions focus on additional, physically plausible conditions. These conditions, which are intended to capture the actual, physical state of affairs, may be imposed on the initial state, the observable, the structure of the system, * rheveling@uos.de † lknipschild@uos.de ‡ jgemmer@uos.de or combinations thereof [14][15][16][17]. In the present paper we primarily discuss results from Ref. [18]. The latter rest on assumptions on all of the above. This paper is organized as follows: In Sect. II we briefly present a main result from Ref. [18]. Furthermore, we elaborate on the lack of predicitive power of this result in cases in which Fermi's Golden Rule applies. In Sect. III we explain some models, each of which comprises a single spin in a magnetic field interacting with a (finite) bath, consisting of spins itself. We initialize the system in a standard system-bath product state and numerically solve the Schrödinger equation, monitoring the system's spin component parallel to the magnetic field. These data unveil the regime of validity of Fermi's Golden Rule with respect to the crucial system parameters. Sect. IV discusses the scaling of the critical interaction strength at which open system predictions start to become unreliable. In Sect. V the implications of the numerical findings from Sect. III on the assumptions and statements from Ref. [18] are named and explained. Eventually, we sum up and conclude in Sect. VI.

II. GARCÍA-PINTOS BOUND AND FERMI'S GOLDEN RULE
To begin with, we state a main result of Ref. [18] (hereafter called the García-Pintos bound (GPB)) in a comprehensive form. The GPB addresses an equilibration time T eq . To further specify T eq we introduce some notation. Let ρ be the initial state of the system. Let furthermore A(t) denote an observable A in the Heisenberg picture and A(t) := Tr[A(t)ρ] its time dependent expectation value. Due to the closed system dynamics being unitary (and the system being finite), A(t) has a well defined "infinite time average" A := A eq , which is routinely considered as the equilibrium value of A in case the observable A equilibrates at all [9]. Consider now a deviation D(t) of the actual expectation value from its equilibrium, i.e., D(t) := ( A(t) − A eq ) 2 /4||A|| 2 , with ||A|| being the largest absolute eigenvalue of A. Consider furthermore an average of D(t) over the time interval [0, T ] denoted by D T . The condition that defines T eq is that D T 1 must hold for T T eq (for nonequilibrating systems such a T eq may not exist [9]). The GPB is an explicit expression for such a T eq (see Eq. (4)), based on ρ, A(0) and H, where H is the Hamiltonian of the system. As the GPB involves somewhat refined functions of the above three operators, we need to specify these before stating the GPB explicitly. A central role takes a kind of probability distribution p jk which is defined as where E j , E k are energy eigenvalues corresponding to energy eigenstates |j , |k . Furthermore, matrix elements are abbreviated as ρ jk := j|ρ|k , A jk := j|A(0)|k . While the GPB is not limited to this case, we focus here on p jk which allow for a description in terms of a probability density function w(G). All examples we present below conform with such a description and it is plausible that this applies to many generic many-body scenarios. Prior to defining w(G), we define w(G, ) as where Θ is the Heaviside function. This is the standard construction of a histogram in which the p jk are sorted according to their respective energy differences E j − E k . It is now assumed that there exists a range of (small but not too small) such that w(G, ) is essentially independent of variations of within this range. The w(G, ) from this "independence regime" are simply abbreviated as w(G). Let the standard deviation of w(G) be denoted by σ G . Let furthermore w max denote the maximum of w(G). The quantities a and Q that eventually enter the GPB are now defined as We are now set to state the GPB: Obviously, the GPB links T eq to the initial "curvature" of the observable dynamics ∂ 2 t A(t) | t=0 (which is practically accessible, cf. Fig. 7). An actual, concrete bound on the equilibration time by means of T eq , however, only arises from Eq. (4) if the numerator can be shown to be in an adequate sense small or at least bounded. This is a pivotal feature on which the "predictive power" of the GPB hinges. The crucial quantities in the numerator are a and Q. As it is practically impossible to calculate a from its definition for many-body quantum systems, García-Pintos et al. instead offer an assumption. They argue that a ∼ 1 may be expected for w(G) that are "unimodal". Unimodal means that w(G) essentially consists of one central elevation like a Gaussian or a box distribution, etc. Indeed, a is invariant with respect to a rescaling as w(G) → sw(sG), as it would result from rescaling the Hamiltonian as H → sH (here s is some real, positive number). García-Pintos et al. also offer various upper bounds on Q for different situations. In the remainder of this section, we explain in which sense the conclusiveness of the GPB is in conflict with Fermi's Golden Rule (FGR). Let us stress that this conflict does not concern the validity or correctness of Eq. (4) as such, the latter is undisputed. It only concerns the assumptions on a and Q, which are required to find an actual value or estimate for T eq . (Note that there is some evidence (cf. Sect. V) that specifically the assumption on a is violated, rather than the assumption on Q). Consider an Hamiltonian consisting of an unperturbed part H 0 and a perturbation H int .
Consider furthermore an observable A, which is conserved under H 0 , i.e. [A, H 0 ] = 0. If H 0 has a sufficiently wide and dense spectrum and λ is small, FGR may apply under well investigated conditions [19][20][21]. The applicability of the FGR approach yields, in the simplest case, a monoexponential decay, i.e.
where τ rel := rλ −2 and r is a real, positive number depending on H 0 and H int . More refined approaches, such as the Weisskopf-Wigner theory or open quantum system approaches, also arrive at such exponential decay dynamics [22,23]. In the relevant case ∂ t A(t) | t=0 = 0, obviously Eq. (6) cannot apply at t = 0. In this case Eq. (6) is meant to apply after a short "Zeno time" τ zeno that is often very short compared to the relaxation time τ rel [21]. (Note, however, that the denominator of Eq. (4) addresses a time below the Zeno time, if the latter is nonzero). We now aim at finding the principal dependence of quantities in Eq. (4) on the interaction strength λ. While the definition of T eq as given at the beginning of the present Sect. does not fix the relation of τ rel and T eq rigorously, for exponential decays it appears plausible to require at least For the denominator of Eq. (4) we find with Eq. (5) for the numerator of Eq. (4). Obviously, the numerator of Eq. (4) diverges in the limit of weak interactions, i.e.
λ → 0. The latter holds even if c 1 = 0. This contradicts the central assumption behind the GPB as outlined below Eq. (4). Hence, the validity of FGR in the weak coupling limit and a conclusive applicability of the GPB are mutually exclusive. This is the first main result of the present paper. Although the practical success of FGR is beyond any doubt, the theoretical applicability of FGR rests on various assumptions on the system in question, so does the applicability of standard open system methods. In order to learn about the applicability of either the GPB or FGR from considering examples, we analyze some spin systems in the following Sect. III by numerically solving the respective Schrödinger equations. This analysis is comparable to numerical investigations performed in Ref. [18]. However, other than García-Pintos et al. we analyze the weak coupling limit and consider system sizes that are too large to allow for numerically exact diagonalization of the respective Hamiltonians.

III. NUMERICAL SPIN-BASED EXPERIMENTS PROBING EQUILIBRATION TIMES
While we analyze a number of concretely specified models below, it is important to note that these models just represent some generic instances of the system-bath scenarios which are routinely considered in open quantum system theory. The (non-integrable) baths share some properties with standard solid state systems, like periodicity and locality (in this respect they differ from the otherwise comparable models addressed in Refs. [24][25][26]). Other than that, the details of our modeling are not peculiar at all. We varied details of the bath Hamiltonians in piecemeal fashion and found all below results unaltered (cf. App. C). Our archetypal model is an isotropic spin-1/2 Heisenberg system consisting of a single system spin coupled to a bath. The bath is rectangularly shaped with 3 × L spins and features periodic boundary conditions in the longitudinal direction resulting in a wheel-like structure (cf. Fig. 1). Thus, the total number of spins is given by N = 3L + 1. The single system spin is subject to an external magnetic field in the z-direction and interacts with three neighboring bath spins in the transverse direction. This model is non-integrable in the sense of the Bethe-Ansatz. The bath Hamiltonian reads where S x,y,z i,r are spin-1/2 operators at site (i, r) and L + 1 ≡ 1. The exchange coupling constant J as well as are set to unity. The Hamiltonian of the system is given by where S x,y,z sys denote the spin-1/2 operators of the additional system spin and B = 0.5. The interaction between bath and system is described by the Hamiltonian and contributes with a factor λ to the total Hamiltonian The considered initial states are product states of a system state π ↑ and a bath state π E,δ . This corresponds to a situation where system and bath are initially uncorrelated and then brought into contact via H int at t = 0. The system state is a projector onto the S z sys -eigenstate corresponding to spin-up. The bath state π E,δ is a projector onto a (small) energy window of width δ centered around a mean energy E.
Concretely, we fix the width of the energy window δ = 0.1, which is very small compared to the scale of the full energy spectrum of the bath. Given the size of the systems it comprises nevertheless a very large number of energy eigenstates. To keep track of finite size effects we increment the baths circumference L in steps of size one, thus adding three spins to the bath in each step. An inverse temperature β is defined as where Ω(E) is the density of states of the bath at energy E. This "microcanonical" definition of temperature is also employed in Ref. [18]. For comparability of different bath sizes we aim at keeping β fixed while incrementing the bath size. As H bath is local, the bath energy is expected to scale linearly with the bath size. Hence, we choose a scaling of the initial bath energy as E ≈ −0.15(N − 1), which corresponds to choosing β ≈ 0.4. Given these specifications of the Hamiltonian and the initial state, we numerically solve the corresponding Schrödinger equation and monitor the expectation value of the z-component of the magnetization of the system-spin, i.e. S z sys (t) . Some results are displayed in Fig. 2 for a schematic overview. In accord with open quantum system theory, these results suggest to distinguish three cases. i. non-Markovian regime: For strong coupling (e.g. λ = 1.0) the magnetization quickly decays to the equilibrium value, i.e.
S z sys mc , in a nonexponential way. The description of these dynamics requires the incorporation of memory effects in some way. While this is a very active field in open quantum theory, we do not investigate this regime any further in the present paper. ii. Markovian regime: For weak coupling (e.g. λ = 0.1) there is a monoexponential decay to the thermal equilibrium.
This exponential decay is in full accord with FGR. A large number of systems ranging from quantum optics to condensed matter fall into this regime [23,27]. It also largely coincides with the field of quantum semi-groups and the Lindblad approach. iii. superweak coupling regime: For very weak coupling (e.g. λ = 0.01) the magnetization does not decay to the thermal equilibrium value at all, it rather gets stuck at a value closer to the initial value, which indicates the breakdown of FGR. This value depends on the interaction strength and on the bath size. In accord with standard open quantum system theory, our below results indicate that this regime only exists for finite baths. We are not aware of any systematic approach to this regime in the literature to date. As the conflict between the GPB and FGR arises in the limit of weak interactions, cf. Eq. (9), we are primarily interested in the transition from the Markovian to the superweak regime. A prime indicator of superweak dynamics is, as mentioned above, the fact that S z sys (t) no longer decays down to the microcanonical expectation value S z sys mc = −0.05 as it does in the non-Markovian and the Markovian regime. Fig. 3 shows the longtime average value of the magnetization plotted over the interaction strength λ for various bath sizes. For sufficiently strong coupling the magnetization decays to the thermal equilibrium value for all bath sizes. For each bath size there exists a critical interaction strength λ crit , below which the magnetization gets stuck at a non-thermal longtime average value, thus signaling the transition from the Markovian to the superweak regime. This critical interaction strength λ crit decreases with bath size. We chose S z sys = −0.04 (horizontal grey line) to define λ crit . It turns out that the below scaling of λ crit is rather insensitive to the exact positioning of this threshold, as long as it is sufficiently close to the thermal equilibrium value. Obviously, one expects S z sys → 0.5 as λ → 0.  Hence, for all mesoscopic to macroscopic systems, and even more so in the thermodynamic limit, the transition to the superweak regime practically never occurs, such that behavior other than Markovian can hardly be expected even for physically very weak interactions. While this finding is another main result of the quantitative analysis at hand, it qualitatively hardly comes as a surprise in a larger context, given the practical success of Markovian quantum master equations. However, to elaborate on this result somewhat further, we present a theory that captures the data in Fig. 4 rather accurately in Sect. IV.
Next we confirm the validity of FGR in the Markovian regime and discuss relaxation/equilibration times in all regimes. The motivation for the latter is twofold: On the one hand equilibration times enter the GPB (cf. Eq. (4)), on the other hand the scaling of equilibration times with the interaction strength may serve as a additional, quantitative indicator for the validity of FGR. Fig. 5 displays the observable dynamics S z sys (t) for 11 randomly selected interaction strengths and bath sizes from the Markovian regime, which is lower bounded by λ crit as determined from Fig. 4 and upper bounded by λ non-Mark ≈ 0.3 for all N . Note the the time axis is scaled with the squared interaction strength such that a collapse of the data onto one decaying exponential indicates the accordance with Eq. (6) and hence FGR. This collapse is evident. In Fig. 6 the relaxation time τ rel is plotted over λ −2 .
Here τ rel is the time at which the magnetization has decayed to 1/e of its original value relative to the equilibrium value (cf. Eq. (6)). In the Markovian regime, i.e. for λ −2 non-Mark ≤ λ −2 ≤ λ −2 crit , the relaxation time scales as τ rel ∼ λ −2 , as predicted by FGR, which also confirms the applicability of FGR in the Markovian regime. At very small λ, i.e. in the superweak coupling regime, τ rel first increases more slowly with increasing λ −2 and likely eventually even decreases to zero. From the data displayed in Fig. 6 this behavior is, however, qualitatively only visible for N = 16 due to numerical limitations at extremely small λ.  Eventually, we directly numerically probe the connection between short time and long time dynamics suggested in Eq. (4). To this end we compute the "initial curvatures" |∂ 2 t A(t) t=0 | for various interaction strengths and systems sizes. The result is displayed in Fig. 7. As expected from Eq. (8), and in full accord with a corresponding statement in Ref. [18], the square root of the initial curvature scales linearly with λ and is practically independent of the system size. We are now set to assess the crucial numerator from Eq. (4) numerically. From Eqs. (4, 7) follows The lower bound to the numerator is displayed in Fig. 8.
Recall that for a conclusive application of the GPB this numerator must be appropriately upper bounded. Correspondingly, Ref. [18] offers estimates for both a and Q. While a ∼ 1 is simply traced back to the unimodality of w, the discussion on the order of magnitude of Q is quite involved. However, in the case of weak interactions, a microcanonical initial bath state comprising a large number of energy eigenstates, and an exponentially growing density of states in the bath (with an exponent β which is not too large), Q may also be expected to be of order unity, according to Ref. [18]. All these conditions apply to the models at hand. However, quite in contrast we find that the numerator grows at least up to values of ca. 55 already for N = 25 and interactions within the range of our numerical accessibility. Moreover, the data do not indicate any "nearby" upper bound of the numerator at 55. This is at odds with a conclusive application of the GPB and another main result of the present paper. Large numerators must occur, as outlined in Sect. II, for large systems in the Markovian regime on the verge to the superweak regime. Although Fig. 8 indicates that in the outer mathematical limit (which may be considered to be physically less relevant) λ → 0, the numerator may eventually be of order unity, its growth appears to continue substantially into the superweak regime. While we are unable to verify this directly, it appears plausible that the unbound growth of the numerator is due to an unbound growth of a, while Q ∼ 1 may very well hold. We argue for this plausibility in Sect. V. We sum up this section as follows: In a potentially wide regime of interaction strengths, which is lower bounded by λ crit , FGR is found to apply as suggested by standard open quantum system theory in the Markovian regime.
The lower bound appears to decrease rapidly with system size N (cf. also Sect. IV), making it practically irrelevant for mesoscopic and macroscopic systems. This relates to the GPR inasmuch as the validity of FGR implies the breakdown of the practical applicability of the GPB at sufficiently weak interactions. We numerically confirmed the occurrence of this breakdown directly for a system comprising N = 25 spins.

IV. GENERAL SCALING OF λ crit WITH TOTAL SYSTEM SIZE
While the results on λ crit in Fig. 4 are model dependent, a similar scaling may be expected whenever H int complies with the eigenstate thermalization hypothesis. This claim is substantiated in the following. The starting point is the assumption that within the superweak regime, the overlap of the eigenstates of the uncoupled system |n, sys + bath and those of the full system |n, sys + bath + int is relatively large, i.e. | n, sys + bath|n, sys + bath + int | ≈ 1 .
A strong indication for this to occur arises from the leading order contributions to a perturbative correction to the eigenstates being small. This condition may, according to textbook level perturbation theory, be approximated as where Ω (N, β) is the density of states of a system comprising N spins (or other similar subsystems) at the energy that corresponds to the inverse temperature β. In Eq. (17) it is assumed that the level spacing within the relevant energy regime may be approximated as being constant. In this case the "mean" level spacing is given by 1/Ω(N, β). Following the eigenstate thermalization hypothesis ansatz [28] we furthermore assume that there exists a typical value for the absolute squares of the matrix elements of the coupling operator, which varies with the energies E n , E m only on energy scales much larger than the one relevant here. Furthermore, the eigenstate thermalization hypothesis suggests a specific scaling of these matrix elements with the density of states. Following the eigenstate thermalization hypothesis we thus assume | m, sys + bath|H int |n, sys + bath | 2 ≈ C 1 Ω (N, β) , where C 1 is some real constant. Exploiting this, Eq. (17) turns into As the sum assumes the finite value π 2 /3, we conclude for the scaling of λ crit , where C 2 is a constant whose concrete value depends on H sys , H bath and H int . Now we turn to an estimate for Ω (N, β). For a sufficiently large Heisenberg spin system (or any other system consisting of N similar, similarly and locally interacting subsystems) it is reasonable to assume that the density of states is Gaussian with mean zero and variance proportional to the particle number N .
Using β = ∂ E log Ω leads to Inserting this into Eq. (20) and fitting for C 2 and α yields the dashed line in Fig. 4, which matches the data quite well. This remarkable agreement in turn backs up the argumentation which lead to Eq. (20). As the density of states is routinely expected to scale exponentially with the size of the system, Eq. (20) indicates that λ crit will generally be exponentially small in the system size and thus the superweak regime will practically never be observed.
V. GARCÍA-PINTOS BOUND AND EXPONENTIALLY DECAYING OBSERVABLES As explained in Sect. II the applicability of the GPB hinges on the two parameters a and Q, both of which should be of order unity to establish a meaningful relation between the short time dynamics and the equilibration time in the sense of the GPB. However, for some of the numerical examples considered in Sect. III at least one of the parameters must be substantially larger than unity. While we are unable to perform a direct numerical check for large system sizes, we strongly suspect that a ∼ 1 is violated at weak interactions even though the corresponding w(G) (cf. Sect. II) is strictly unimodal. In the remainder of the present section we explain and back up this claim. Consider the mathematically simple case of an infinite temperature environment in the initial state ρ = (S z sys + 1 sys /2) ⊗ 1 bath /d bath . This initial state yields S z sys (t) = Tr{S z sys (t)ρ} = Tr{S z sys (t)S z sys } for the dynamics of the observable S z sys (t) , which may be rewritten as Now consider the distribution p jk as defined in Eq. (1) for this initial state and choice of observable, i.e., A = S z sys .
Due to the system being non-integrable in the sense of a Bethe ansatz, the eigenstate thermalization hypothesis may be expected to hold, yielding | j|S z sys |j | 2 ≈ 0. Exploiting this, the insertion of Eq. (25) into Eq. (24) yields To the extend to which p jk may indeed be replaced by a smooth probability density as discussed around Eq. (2), Eq. (26) may be rewritten as Thus, for the present scenario, w(G) is essentially the Fourier transform of the observable dynamics S z sys (t) . Based on the numerical findings displayed, e.g., in Fig.  2 it appears plausible that S z sys (t) will be an exponential decay for infinite temperature initial states as well. Therefore, w(G) will be Lorentzian. While a Lorentzian distribution is clearly unimodal with one well-behaved maximum, its variance diverges. Consequently a, as defined in Eq. (3), diverges as well. Thus, in contrast to the assumptions in Ref. [18], a ∼ 1 does not hold. This is the last main result of the present paper. Of course σ G cannot really diverge in any system featuring a finite energy spectrum. However, the finiteness of the spectrum essentially causes a cut-off of the tails of the Lorentzian at some frequency. This cut-off actually renders the standard deviation σ G finite. Nonetheless, this standard deviation does not reasonably reflect the width of w(G). It will be much larger than other measures of the width such as the full-width-at-half-maximum, etc. Some attention should also be paid to the question whether or not the GPB scales with the size of the environment. (Earlier works presented upper bounds that explicitly depend on the size of the environment, which is often seen as a drawback [10,11]). While the GPB does not explicitly dependent on the size of the environment, the latter may enter via the parameter a. For any (weak) interaction strength λ there exists a system size N (λ crit ) above which FGR applies. Above that size the GPB is independent of N . Below or at N (λ crit ), however, the numerator in Eq. (4) may depend on N rather strongly. For arbitrarily small λ this N (λ crit ) may become arbitrarily large. Thus, in the class of models discussed in the paper at hand, one can always find instances for which the GPB depends on system size even for very large systems, i.e. N 1.

VI. SUMMARY AND CONCLUSION
In the paper at hand we conceptually and numerically analyzed an upper bound on equilibration times presented in a recent paper by García-Pintos et al. To this end, we investigated a standard system-bath setup by monitoring the system's magnetization for various bath sizes and interaction strengths. This numerical investigation is based on the solution of the time dependent Schrödinger equation for the full system, including the bath. We identified a Markovian regime of interaction strengths λ in which Fermi's Golden Rule holds, i.e., the system thermalizes in an exponential way and the equilibration time scales as λ −2 . This relates to the García-Pintos bound inasmuch as the validity of Fermi's Golden Rule and the usefulness of the García-Pintos bound are analytically shown to be mutually exclusive at sufficiently small λ. At extremely small λ, we indeed find a "superweak" regime in which Fermi's Golden rule does not apply. This regime (in principle) exists for finite baths and is reached below some λ crit which is shown to scale inversely exponentially in the bath size, suggesting that the superweak regime practically ceases to exist when considering moderately large systems. However, in the superweak regime the García-Pintos bound may eventullay regain applicability, although its non-applicability is found to extend also into the superweak regime.

ACKNOWLEDGMENTS
Stimulating discussions with R. Steinigeweg and J. Richter are gratefully acknowledged. We also thank L. P. García-Pintos and A. M. Alhambra for interesting discussions, as well as A. J. Short for a comment on an early version of this paper. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) within the Research Unit FOR 2692 under Grant No. 397107022.
Details of our numerical implementation will be discussed in this appendix. We make use of the concept of typicality (cf. App. A) and use a time evolution algorithm (real and imaginary) based on Chebyshev polynomials (cf. App. B). Since the full Hamiltonian conserves magnetization, we perform the procedure outlined below in each magnetization subspace. The dynamics in the full Hilbert space are obtained by piecing together the contributions of each magnetization subspace weighted with their respective binomial weight. Note that the main hindrance to our calculations is not the exponentially large Hilbert space dimension, but rather the extremely long times that have to be reached in real time. In App. C a result for randomized bath couplings is shown.

A. Typicality
As mentioned in Sect. III, the initial state is a product state of a microcanonical bath state and a projected spin-up system state. Since the numerical integration of the von-Neumann equation can be cumbersome, we make use of the concept of typicality, which states that a single "typical" pure state can have the same thermodynamic properties as the full statistical ensemble [29]. Not only is it more memory efficient to work with pure states, the availability of efficient time evolution algorithms for pure states, e.g. Runge-Kutta or Chebyshev polynomials, constitutes a major advantage. To find such a typical state, a pure state |φ is drawn at random from the Hilbert space according to the unitary invariant Haar measure. It can be shown [30] that for the overwhelmingly majority of random states |φ , the pure state |ψ exhibits effectively the the same thermodynamic behavior as the mixed state ρ, i.e.
Importantly, the induced error = (|ψ ) has mean zero, i.e. = 0, and a standard deviation that scales inversely proportional to the square root of the effective Hilbert space dimension, i.e. σ( ) ∝ 1/ √ d eff [31]. The effective dimension d eff = 1/Tr{ρ 2 } is a measure of how many pure states contribute to the mixture ρ.
In the paper at hand the initial state ρ (cf. Eq. (14)) is a projection operator. Therefore, it is permissible to drop the square root in the numerator in Eq. (A.2). Now the projectors π ↑ ⊗ 1 and 1 ⊗ π E,δ need to be applied to the state |φ , which is an element of the product Hilbert space. As we are working in the Ising basis, the action of π ↑ ⊗ 1 on |φ is easily implemented by setting corresponding components of the state vector to zero. Since it is unfeasible to diagonalize the full many-body Hamiltonian, we replace the bath projector by a Gaussian filter which suppresses contributions of energy eigenstates far away from the desired energy E, resulting in a narrowly populated energy window of width (variance) δ.
The Gaussian filter is applied with a Chebyshev algorithm (cf. App. B). Lastly, the resulting wave function is normalized to obtain the state |ψ as in Eq. (A.2). In this scenario the effective dimension d eff is essentially the number of states in the energy window. Increasing the size of the system, while keeping δ fixed, results in an exponentially growing effective dimension d eff and therefore in a negligible typicality error for moderately sized systems.

B. Chebyshev polynomials
A Chebyshev type algorithm is employed in order to evolve a pure state in real and imaginary time [32,33]. Say it is desirable to approximate a scalar function f (x) in the interval [−1, 1] by a polynomial expansion, i.e.
f (x) ≈ n c n P n (x) , (B.1) with coefficients c n and polynomials P n of order n. The unique set of polynomials that minimizes the maximum error in this interval is called Chebyshev polynomials of the first kind. They are denoted by T n (x) and can be written down recursively as with T 0 (x) = 1 and T 1 (x) = x. They are orthogonal with respect to the weighted scalar product with C 0 = 1 and C n>0 = 1/2. In order to approximate the time evolution operator, the bandwidth of the Hamiltonian has to be rescaled accordingly. Defining a = (E max − E min )/2 and b = (E max + E min )/2, where E max (E min ) is the maximal (minimal) energy eigenvalue, the rescaled Hamiltonian is obtained byH = (H − b)/a.
In practice a small safety parameter is chosen that ensures that the rescaled spectrum lies well within While the investigated model class may seem peculiar, it has just been chosen as one generic representative of the whole of condensed matter type systems. To exclude that overall results are just due to any unintentional, subtle conserved quantities, etc., we redid parts of our numerical analysis with randomized bath couplings, drawn from a Gaussian distribution with mean zero and standard deviation 0.2. One result is displayed in Fig. 9 for N = 25 and λ = 0.2. One readily verifies that the curves coincide nicely. This hints at the generic nature of our model class. Comparison between setups with couplings set to unity and randomly drawn couplings. There is no apparent difference.