Space-local memory in generalized master equations: Reaching the thermodynamic limit for the cost of a small lattice simulation

Lattice models play a key role in understanding physical and chemical phenomena. For instance, the Holstein^1,2 and Fröhlich³ models shed light on polaron formation and electrical transport in semiconductors,^4,5 the Hubbard model⁶ helps elucidate the mechanisms of high-temperature superconductivity,⁷ and the Ising model⁸ is used to interrogate magnetism⁹ and phase transitions.¹⁰ However, while modern algorithms can efficiently simulate the quantum dynamics of small lattices over short times,^11–28 reaching sufficiently large systems and long times to compare with experiments remains a fundamental challenge. This is because these methods often scale exponentially or, at best, polynomially with lattice size and simulation time, rendering the thermodynamic limit inaccessible.

The severity of this limitation becomes clear when calculating dynamic properties, e.g., conductivities, viscosities, and diffusion constants, which are sensitive to finite-size effects.^29–36 For example, finite-size effects can cause simulations to underestimate the diffusion constants of polymers near the glass-transition,³⁷ relaxation times of glass-forming liquids,³⁸ the Curie temperature for Ni nanoparticles,³⁹ and the diffusion constant and viscosity of model fluids;⁴⁰ overestimate the critical fermion–phonon coupling, causing the metal-to-Peierls phase transition in a Holstein–Hubbard lattice;⁴¹ and yield apparently non-converging mobilities of dispersive Holstein polarons.⁴² These examples reveal the need for computing the dynamics of lattice models in thermodynamically large systems over long timescales.

Generalized Master Equations (GMEs) have emerged as a powerful tool for reducing the computational cost of dynamical simulations.^43–62 GMEs are exact non-Markovian equations of motion for nonequilibrium averages, correlation functions, and even multi-time correlators of select variables that encapsulate the effects of an environment into a memory kernel.^63–65 In dissipative systems, the memory kernel decays to zero over a finite memory lifetime, which can be shorter than the relaxation time of the desired correlation function. Thus, in principle, one can use a reference simulation over the memory lifetime to construct a GME that predicts the dynamics of the desired correlation function to arbitrarily long times. This temporal truncation of memory at its lifetime can reduce the computational cost of simulating the quantum or classical dynamics of complex many-body systems in different problems, including charge transfer reactions in solution,^66,67 protein folding,^68–70 nonlinear spectroscopy,^61,71 and transport.^72,73 However, to construct a GME from a short-time reference simulation, it must satisfy two conditions: (1) the simulation time must span the memory kernel lifetime and (2) the reference calculation must be performed in the same system whose dynamics one intends to interrogate with the GME. If one constructs a GME using a small lattice simulation, particles encounter the lattice boundaries and manifest finite-size effects: one reduces the cost but still obtains the wrong answer. Hence, one must be able to afford an admittedly short-time reference simulation, but of a thermodynamically large lattice. The poor scaling of dynamical methods with system size renders this calculation at best impractical and at worst impossible.

Here, we propose a novel approach to lattice problems that exploits our observation that certain GME formulations display a finite spatial memory to motivate truncating memory in both time and space. This allows us to employ short-time reference simulations of small lattices to generate the exact quantum dynamics of thermodynamically large lattices over arbitrarily long times, at the cost of only a small reference calculation. We demonstrate the strengths of this method by applying it to nonequilibrium polaron formation and transport in dispersive Holstein lattices. Enabled by our space-local GME, we simulate, for the first time, the exact nonequilibrium quantum dynamics of small polaron formation, relaxation, and flow in thermodynamically large one-dimensional (1D) and two-dimensional (2D) lattices with up to 900 sites over 100 ps, free of finite-size contamination. Our method is model-agnostic and can be expected to enable the efficient investigation of nonequilibrium excitation dynamics in dissipative lattices displaying local interactions.

The dispersive Holstein model offers a physically transparent description of small polaron formation and transport in many semiconductors that exhibit strong carrier–phonon interactions, including organic crystals,⁷⁴ polymers,⁷⁵ and nanomaterials.⁷⁶ This model describes the migration of charge carriers (electrons or holes) or excitons, which interact with their environment, influencing local lattice (nuclear) motions. This interaction causes the material to deform around the carrier, forming a polaron.

Consistent with previous works,^5,42,72 we focus on the dilute limit (one-carrier manifold) of homogeneous lattices, making all our parameters site-independent. Without loss of generality, we set ϵ_i = 0, the nearest-neighbor hopping integral v_ij = v, inter-site distances r₀ = 5 Å, and all spectral densities to an Ohmic–Debye form commonly used to mimic dissipation in the condensed phase,⁸²  $J(ω)=ηωcω/(ω2+ωc2)$⁠. Here, η/2 is the lattice reorganization energy of an occupied lattice site and 1/ω_c is the phonon environment’s decorrelation time.

In our simulations, we track the population matrix C(t), where $Ci,j(t)=Tr[âi†âie−iLtρ̂j(0)]$ with the initial condition $ρĵ(0)=âj†âje−βĤph/Tr[e−βĤph]$⁠, which corresponds to creating an excitation at site j at t = 0 when all the phonon modes are in thermal equilibrium with the electronic ground state at inverse temperature β = 1/k_BT. C_i,j(t) quantifies the probability of finding the polaron at site i at time t given that it was initiated at site j. In homogeneous systems with equivalent sites and periodic boundary conditions, C(t) exhibits translational invariance, i.e., C_i,i+lmod N(t) = C_j,j+lmod N(t) for all i, j, l. This allows one to identify all distinct elements of this matrix using their relative site index, |i − j| = k, allowing one to rewrite C_i,j(t) = C_|i−j|(t) = C_k(t). Hence, one only needs to perform a single simulation starting the carrier at any particular site, j, and record the time-dependent probability of finding the carrier at all lattice sites, i, as a function of time to construct C(t).

We begin with the $dMSDdt$ dynamics of dispersive Holstein polarons on a 1D lattice. To align with previous work on organic semiconductors,^42,86 we take η = 323 cm⁻¹, v = 50 cm⁻¹, and ω_c = 41 cm⁻¹. Beyond encoding the diffusion constant, $dMSDdt$ provides insights into polaron formation, the transition from nonequilibrium relaxation to diffusive transport, and the onset of finite-size contamination.⁴² For example, in the reference (HEOM) dynamics for the 20-site lattice (solid pink line) in Fig. 1(a), the initial hump at ∼100 fs indicates far-from-equilibrium lattice reorganization due to polaron formation. This is followed by a plateau that suggests diffusive motion between 800 and 2000 fs. While this region is not truly flat, it approximates the diffusion constant.⁴² This putative plateau falls off at ∼2000 fs, marking the onset of finite-size effects. We denote this onset timescale as τ_R.

We address this fundamental problem by leveraging the concept of spatial memory to develop a space-local (SL) GME that uses only a short-time simulation of a small lattice to predict the dynamics of thermodynamically large lattices over arbitrary times. We begin by noting that, like C(t), $U(t)$ exhibits translational invariance in homogeneous lattices, allowing us to use relative indices to identify the elements. Figure 2(a) shows that as inter-site distance increases, the elements of $U(t)$ become progressively smaller, $Ud+2<Ud+1<Ud$⁠, suggesting that our TL-GME generator exhibits decaying spatial memory. Furthermore, one may posit that generator elements connecting sites separated by a distance greater than a characteristic memory distance, $dU$⁠, become negligibly small, allowing one to set them to zero. We illustrate this concept of spatial memory truncation in Fig. 2(b), where we discard the elements of $U(t)$ beyond a particular distance, $d=dU$⁠, as they become negligibly small.

We turn to testing the validity of SL-GME. We begin by examining the ‖L‖₂ error metric for the SL-GME as a function of the proposed memory distance cutoff, d, to identify the characteristic memory distance, $dU$⁠. Figure 3(a) shows that the SL-GME reproduces the reference dynamics (solid pink line) when $dU=3r0$ (dotted brown line), but not for smaller cutoff distances. The plateau in the error as a function of the distance cutoff in Fig. 3(b) confirms that $dU=3r0$⁠, achieving an average error of less than 1%. These results thus support the proposal of a memory distance cutoff that preserves accuracy in the SL-GME dynamics.

Since obtaining the generator up to its lifetime $τU$ allows one to simulate the dynamics for all time, one may hypothesize that constructing the generator up to its characteristic memory distance $dU$ should enable one to simulate a lattice of any size, albeit of the same dimension. This is the central insight of our space- and time-local GME (STL-GME). To achieve this in a 1D lattice, we propose augmenting the dimension of $U(t)$⁠, which is a 3-tensor with dimensions (N, N, tsteps), to (M, M, tsteps), where M > N and all new elements are assigned a value of 0. Furthermore, truncating the temporal dimension of our augmented generator at time $τU$ enables us to easily generate the dynamics after $τU$ via simple matrix multiplication. This spatial and temporal truncation of $U(t)$ followed by its spatial augmentation constitutes our STL-GME. We offer details of its implementation in  Appendix C. The resulting STL-GME should thus offer a route to employ the dynamics of a small lattice over short times to simulate the behavior of a thermodynamically large lattice across arbitrarily long timescales.

We test our STL-GME’s performance by turning again to the Holstein lattice in Figs. 1 and 3. Since $U$ is an intrinsic property of the dynamics, one can extract it from short-time $t≤τU$ reference dynamics of a lattice with $N≥2dU+1$ sites. One can then employ $U(t)$ to predict relaxation dynamics over arbitrarily large spaces and times, removing all finite-size artifacts. We validate this idea in Fig. 4, where we construct the generator for our STL-GME dynamics for a lattice with N = 20 sites using a reference simulation of a lattice with only $N=8>2×(dU=3)+1$ sites (solid purple line). Our size-augmented STL-GME in Fig. 4 (light purple dots) reproduces the reference dynamics for a lattice with N = 20 sites (solid pink line). In fact, the STL-GME can access the dynamics of arbitrarily large systems, indefinitely delaying the onset of finite-size effects. As we show in the inset of Figs. 4, a 20-site lattice exhibits finite-size artifacts around ∼2000 fs, whereas our STL-GME calculation for the 100-site system shows no signs of finite-size effects, even by 25 ps.

Yet, the fact that N can be so small in the reference calculation used to construct the generator in Fig. 4 is surprising. Here, $τU≃800$ fs, but for an 8-site lattice, finite-size effects manifest in dMSD/dt at ∼300 fs, i.e., τ_R ≃ 300 fs. However, despite $τU>τR$⁠, the STL-GME and reference calculations agree, meaning that the STL-GME correctly simulates the larger lattice and is not contaminated by the smallness of the reference system.

How can one understand this seeming contradiction where the generator lifetime, $τU$⁠, can be longer than τ_R, the onset of finite-size effects in dMSD/dt? Immediately, at τ_R, the finite-size effect contaminates only the elements of $U$ connecting the origin to the most distant sites. This is true because the population density is conserved and its current local. In time, this artifact spreads to elements of $U$ progressively closer to the origin, as the density interferes with itself. One might hypothesize that finite-size effects only manifest in the STL-GME when the elements of $U$ that survive the spatial truncation have been contaminated—a time longer than τ_R.

To confirm this hypothesis, in Fig. 5, we compare the STL-GME generator and dynamics constructed from reference data for two different lattice sizes with varying cutoff distances. The parameter regime is the same as in Fig. 4, where we found $dU=3r0$ and $τR=300$ fs, thus requiring a minimal lattice of N = 8. Keeping $τU=800$ fs constant, we consider two complementary cases. First, for a too-small 6-site lattice, we know that obtaining a result free of finite-size effects is impossible because the periodic boundary is always closer than 3 sites from the initial polaron position. Indeed, $UN=6(d=3r0,t)$ differs significantly from the benchmark $UN=20(d=3r0,t)$ [Fig. 5(a)], with the former producing STL-GME dynamics that deviate from the reference simulation for a 20-site lattice [Fig. 5(b)]. Second, we consider the previously sufficient 8-site lattice but increase the cutoff to $dU=4r0$⁠. Since $N<2dU+1$⁠, we predict finite-size contamination in $UN=8(d=4r0,t)$⁠. This is evident when comparing with the benchmark $UN=20(d=4r0,t)$ [Fig. 5(c)] and from the STL-GME dynamics, which again deviate from the reference [Fig. 5(d)]. That is, while the 8-site lattice generator’s first four entries are correct for this $τU$⁠, the fifth element is contaminated by finite-size artifacts, degrading the predicted dynamics. In both of these pathological cases, truncation in time and space happens after knowledge of the lattice’s finite size reaches those elements in $U$ that we keep. Therefore, a lattice with N = 8 sites converges to the right value when one truncates spatially at $dU=3r0$⁠, as we confirm by comparing $UN=8(d=3r0,t)$ and $UN=20(d=3r0,t)$⁠, and the resulting STL-GME prediction with the reference results for a 20-site lattice, shown in Fig. 4. Hence, the STL-GME is even more efficient in accessing the dynamics of thermodynamically large systems from small reference calculations than one might have expected based on simple physical arguments.

One may further ask if our protocol offers any advantage over describing this transport with a simple diffusion equation. The answer is yes. For example, we have recently demonstrated that, in the case of this small polaron-forming system, nonequilibrium lattice relaxation can exponentially delay the onset of diffusion.⁸⁷ This means that the diffusion equation only becomes valid at exponentially long times, requiring large and expensive calculations, which our STL-GME reduces dramatically. One may, however, argue that the putative plateau region in Fig. 4 around t = 1000–2000 fs provides a sufficiently accurate estimate for the diffusion constant. Obtaining this approximate result requires the simulation of a 20-site lattice over ∼2000 fs. However, recovering the exact value of the diffusion constant from our STL-GME requires only the simulation of an 8-site lattice over ∼820 fs, which reduces the cost by a factor of ∼8. In short, our STL-GME seamlessly bridges numerically rigorous, microscopic quantum dynamics and mesoscale descriptions of transport, like the diffusion equation.

Finally, while we employ the population-based projector, the workflow outlined here is also compatible with other spatially local projectors, like the Argyres–Kelley⁸⁸ one, which is useful when probing decoherence, coherence transfer, and time-resolved spectroscopies. In such cases, one would need N/2 + 1 reference calculations—rather than only one—to construct the GME and its extension. Despite this higher cost, it may be expected to lead to even shorter memory lifetimes.^52,55,56

Thus far, we have assumed only nearest-neighbor coupling, raising questions about the extent to which our conclusions depend on the locality of the Hamiltonian. For example, if the electronic coupling extends further, over second and third nearest neighbors, do the cutoffs extend by 2r₀, or is there a more pernicious, qualitative change?

To probe this potential sensitivity to non-locality, we recalculate our dispersive Holstein dynamics with beyond-nearest-neighbor electronic couplings for the parameter regime of Fig. 3, where we had identified $dU=3r0$ when considering only nearest-neighbor coupling. Motivated by atomic orbital decays, we consider interactions that decrease exponentially with inter-site distance up to the third nearest neighbors, as shown in Fig. 6, with $v(d)=ve−a(d−r0)$⁠, $a=2r0−1$⁠. Even in this non-locally interacting problem, Fig. 6 shows that our STL-GME captures the reference dynamics with $dU=4r0$⁠. This change in cutoff distance of one lattice spacing is even less severe than the 2r₀ one might have anticipated. Therefore, our STL-GME efficiently handles Hamiltonians with non-local, albeit short-range, interactions. While more general models include long-range, through-space effects like carrier–carrier interaction, screening in the condensed phase^89–96 usually allows for a truncation that can be modeled similarly to the modification tested here and should therefore preserve the space-local arguments central to our STL-GME method.

The major advantage of our STL-GME is the significant computational savings it offers. In particular, exact numerical methods scale polynomially or exponentially with N and simulation time. In contrast, our STL-GME scales only as N² with a prefactor that heavily suppresses the cost to keep it effectively constant (see Fig. 7) and (sub)linearly⁹⁷ in time. The main source of computational expense in our STL-GME arises from the short-time reference calculation on a small lattice needed to construct the generator. Figure 7 compares the computational costs of HEOM and the STL-GME built from a small reference calculation as a function of N. We fit the time costs for HEOM with a polynomial regression (gray dashed line in Fig. 7). Even with our highly optimized implementation of HEOM that exploits recent advances in dynamic filtering of auxiliary density matrices⁸³ and the n-particle approximation,⁸⁵ the direct simulation for a 100-site lattice conducted over 25 ps is ∼300 times more computationally expensive than our STL-GME for the same parameter regime. We can further extend the 100-site lattice’s dynamics to arbitrarily long times using our STL-GME for the trivial cost of repeated small matrix multiplications. This STL-GME-enabled simulation allows us to find that finite-size effects emerge at ∼65 ps for the 100-site lattice—a timescale that is currently unreachable for a system of this size with a direct HEOM simulation. Such a calculation would require at least 10 500 CPU hours (optimistically assuming a linear relationship between simulation time and CPU time), making it 750 times more computationally costly than our STL-GME. Instead, with our STL-GME, the same 100-site simulation requires effectively the same amount of time as an 820 fs-long simulation for an 8-site lattice (13.5 CPU hours), as we can compute the STL-GME dynamics in less than a minute.

Our STL-GME method can easily generalize to higher-dimensional lattices. This enables us to access, for the first time, the exact quantum dynamics of polaron formation and transport in dimensions above one. Accessing transport in 2D is particularly significant as it yields the spatiotemporal spread of polarons on a surface—a phenomenon that recent microscopy experiments measure.^98–105 Such simulations can elucidate the microscopic factors that influence energy and charge flow in materials and offer insights on how to modify them for specific applications, such as energy storage or battery development.^106–109 

For our demonstration, we select Hamiltonian parameters characterized by high reorganization energy relative to the hopping integral and a fast phonon decorrelation speed, appropriate for polarons in organic crystals.¹¹⁰ To interrogate the effect of various energy scales on polaron transport, we examine a homogeneous 2D lattice characterized by two types of hopping integrals that can cause anisotropic flow, $V⊥$ and $V‖$⁠, as illustrated in Fig. 8(a). Figures 8(b)–8(d) illustrate the polaron motion on a 30 × 30-site 2D lattice with a total of 900 sites. Such a large-scale simulation is made possible exclusively through our STL-GME method (see  Appendix D for implementation details).

We turn to the effect of varying $V⊥$ relative to $V‖$⁠. Unsurprisingly, when $V⊥=V‖$⁠, we observe a symmetric spread of the polaron at t = 10 ps, as shown in Fig. 8(b). Anisotropic motion can be expected to emerge when $V⊥≠V‖$⁠. Indeed, when $V⊥=2V‖$ and $V⊥=4V‖$⁠, the polaron density spreads more rapidly along the vertical axis, as shown in Figs. 8(c) and 8(d), respectively, with the ratio $V⊥/V‖$ dictating the anisotropy of the polaron distribution on the 2D surface. This demonstrates that STL-GME can efficiently simulate nonequilibrium polaron motion in higher-dimensional systems, even in the thermodynamic limit. Hence, this framework offers an accurate and efficient means to simulate polaron transport over the length and time scales needed to compare to experiments.

Our space locality arguments are also compatible with the time-nonlocal description of memory.^63–65 In  Appendix E, we formulate the space-local time-nonlocal (SL-TNL) GME, where we further unify the SL framework with the transfer tensor method.¹¹¹ Interestingly, we find that the time-local characteristic memory distance, $dU$⁠, is shorter than its time-nonlocal counterpart, $dK$⁠. We also find that the generator and memory kernel lifetimes in the time-local and time-nonlocal formulations follow a similar inequality, $τU≤τK$⁠, consistent with our previous findings in biomolecular dynamics.⁶⁹ 

For systems involving coupling between spins or fermions with phonons, increasing the reorganization energy leads to high electronic–nuclear entanglement. This generally translates into greater computational cost for many numerically exact^11–14,16 and approximate dynamics methods.^112–114 For instance, the number of differential equations that must be propagated at each step in HEOM grows significantly with increasing reorganization energy, raising its computational cost. Similarly, in matrix product state-based propagation schemes, the bond dimension increases with increasing reorganization energy and simulation time,^22,26 further escalating computational demands. In contrast, in our SL-GME approach, increasing reorganization energy reduces $dU$ (see Fig. 9). This allows us to perform reference calculations on smaller lattices, providing dramatic computational savings in systems with large reorganization energies.

While it is well-known that the temporal truncation of memory in traditional GMEs can significantly lower the computational cost of dynamic simulations over arbitrary times, here we introduce space-local GMEs where spatial truncation of memory can dramatically lower the cost of accessing the dynamics of complex many-body systems of arbitrary sizes. By integrating both spatial and temporal truncation in the same GME within a time-local framework, we have developed a novel STL-GME that allows us to leverage the dynamics of small lattices over short times to exactly simulate the dynamics of thermodynamically large systems across arbitrary times, free from finite-size effects. In addition to our time-local formulation, we have developed analogous continuous and discrete time-nonlocal GMEs.

We have demonstrated the benefits offered by our STL-GME when simulating polaron formation and transport in dispersive Holstein lattices in 1D and 2D with short-range couplings. In addition to homogeneous lattices, one can apply our STL-GME to periodic ones.⁸⁷ Our approach, however, can be expected to find applications across a wide range of dissipative lattice problems to interrogate the nonequilibrium relaxation dynamics of charge, electronic energy, and heat flow. We have shown, for example, that our SL-GMEs are broadly compatible with systems displaying short-range interactions and in arbitrary dimensions. In addition, like previous GMEs, our SL-GMEs are agnostic to the choice of dynamics solver (classical or quantum, exact or approximate) used to construct the time-local generator or memory kernel, broadening its applicability to a wide variety of systems and phenomena. While here we use the nonequilibrium population correlation function for the formulation of our STL-GME, our method can also be expected to be compatible with equilibrium correlation functions as long as the operators being correlated are local in space. By offering access to experimentally relevant relaxation processes in thermodynamically large systems over arbitrarily long times, our STL-GME can be expected to provide a transformative tool to uncover and explain new dynamical phenomena in periodic materials.

The donors of the American Chemical Society Petroleum Research Fund are acknowledged for partial support of this research (Grant No. PRF 66836-DNI6). A.M.C. acknowledges the support from a David and Lucile Packard Fellowship for Science and Engineering. S.B. acknowledges the John Bailar Memorial Endowment and the Marion L. Sharrah Fellowship for partial support of the research. T.S. is the recipient of an Early Career Fellowship from the Leverhulme Trust. We thank Professor Qiang Shi for sharing his HEOM code with us. This work utilized the Alpine high-performance computing resource at the University of Colorado Boulder. Alpine is jointly funded by the University of Colorado Boulder, the University of Colorado Anschutz, Colorado State University, and the National Science Foundation (Award No. 2201538).

The authors have no conflicts to disclose.

Srijan Bhattacharyya: Formal analysis (lead); Investigation (lead); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Thomas Sayer: Formal analysis (supporting); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Andrés Montoya-Castillo: Conceptualization (lead); Funding acquisition (lead); Methodology (equal); Project administration (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

HEOM¹¹ is a numerically exact method that predicts the dynamics of the reduced density matrix of an open quantum system by mapping environmental degrees of freedom to auxiliary density matrices (ADMs) that quantify the number of coupled differential equations being solved simultaneously. The reduced density matrix provides information about the polaron population.

We converge all HEOM calculations with respect to the hierarchical depth L, the number of Matsubara frequencies K, and the time step δt. For our 1D simulations, L = 22, K = 1, and δt = 0.25 fs. When changing reorganization energy for Fig. 9, we use L = 26, K = 2, and δt = 0.25 fs. Due to the extreme computational cost, we employed a high-temperature approximation for 2D simulations, setting K = 0, and converged the hierarchical depth and time step to L = 26 and δt = 0.25 fs. In all our simulations, we apply dynamic filtering,⁸³ using a threshold of N_cut = 10⁻⁷ atomic units, in line with previous work.^42,72 We also employ an n-particle approximation⁸⁵ to minimize numerical costs without affecting the accuracy of the results.

Here, we outline our protocol for building a size-augmented STL-GME for a 1D lattice using only a short-time simulation of a small lattice. Here, we use HEOM to generate the reference dynamics, but our method is agnostic to the choice of solver.

1. Construct the generator $U(t)$ from the reference C(t) using Eq. (8).

2. Identify the lifetime $τU$ by computing the ‖L‖₂ error between reference and GME-generated dynamics as a function of the proposed lifetime, τ, $Error(τ)=1Nt‖Cref(t)−CGME(t;τ)‖2$⁠, and identify when the error function plateaus to a sufficiently low value. N_t represents the number of time points predicted by TL-GME, i.e., the number of time points used to construct the generator subtracted from the total number of time points in the reference dynamics. Our threshold per element of C(t) is 3 × 10⁻⁸.

3. Identify characteristic memory distance $dU$ in a similar way by analyzing the error between reference dynamics and GME dynamics via the error as a function of the proposed cutoff distance, d, $Error(d)=1Nt‖Cref(t)−CGME(t;d)‖2$⁠. For SL-GME, N_t is the total number of points. Our threshold per element of C(t) is 6 × 10⁻⁸.

4. If one finds $dU<N/2$⁠, this means the reference simulation is sufficiently large for system size extension via generator augmentation. Keep the generator up to its lifetime $τU$ only. Truncate the generator by setting its entries connecting sites with a relative distance $d>dU$ to 0.

5. Augment the generator’s dimension from (N, N, tsteps) to (M, M, tsteps), where M is the length of the extended lattice. Then, populate all additional new entries in the matrix with 0.

6. Employ $U(t)[M,M]$ to propagate the dynamics. For time $t≤τU$⁠, use Eq. (6). For $t>τU$⁠, use Eq. (7).

In our STL-GME method, we discard the elements of $U$ with relative lattice spacings larger than $dU$ by setting them to 0. Our generator truncation—like singular value truncations in tensor network-based methods¹¹⁷—can cause violations of population conservation. We track this population loss via σ(t) = |1 − ∑_kC_k(t)|. Ideally, σ(t) should be as close to zero as possible. Figure 10(a)—which corresponds to the parameter regime of Fig. 4 where we employ an 8-site reference simulation to predict the dynamics of a 20-site lattice—reveals that our spatial truncation causes a 1% population loss over the first 1 ps, which will continue to grow with simulation time.

To conserve the population, we propose the following two schemes:

1. Redistribution: Add all discarded elements and distribute them equally among all remaining elements of $U$ after truncation,

   $U[i,j;d<dU]=U[i,j;d<dU]+1NdU∑jU[i,j;d>dU],∀i,$

   (C1)

   where $NdU$ is the number of elements in $U$ that we keep as nonzero entries in our generator after spatial truncation.
2. Renormalization: After spatial truncation, renormalize the generator $U$ matrix at each time step,

   $U[i,j;d<dU]=U[i,j;d<dU]∑jU[i,j;d<dU],∀i.$

   (C2)

   Figure 10 shows that when we add our population conservation schemes by redistributing or renormalizing spatially truncating generator $U$⁠, we can maintain the population loss in the order of 10⁻¹¹ and 10⁻¹², even for longer simulation time. Figures 10(b) and 10(c) show the population loss quantity for redistribution and renormalization schemes, respectively. In our study, we employ the redistribution scheme for all 1D simulations and the renormalization scheme for 2D simulations. While we show that we modify the generator to conserve the population, one could, in principle, imagine some level of modification in the population level to conserve the population.

Figure 2(a) shows that increasing lattice distance d makes the elements of the generator $U$ smaller. At some point, the size of these generator elements is commensurate with the statistical error or numerical precision error of the reference dynamical solver. While keeping the full generator intact ensures that the generator recovers the reference dynamics, truncating generator elements that may be at the level of the solver error can introduce small disagreements between the STL-GME and reference dynamics because of imperfect cancellation of error.

To see how this finite precision error in the underlying solver can manifest in the STL-GME, consider generator elements connecting sites separated by inter-site distances of 6r₀ and 7r₀ in Fig. 12 (purple lines), which have contributions on the order of 10⁻⁶–10⁻⁷. In our reference HEOM simulations, from which we build the generator, we use dynamic filtering with a threshold of N_cut = 10⁻⁷ atomic units. Including these small noisy elements in the truncated $U$ risks introducing noise into our dynamical quantities, with derivative quantities magnifying the effect of such noise. In particular, when we choose $dU=6r0$ and compare the exact and STL-GME-predicted MSD, as shown in Fig. 11(a), one cannot notice any visual difference. In contrast, when we compare the more sensitive dMSD/dt, we find a minor noisy behavior in the green-encircled region of Fig. 11(b). This noise arises because elements in $U$ for d ≥ 6r₀ are discarded and set to 0, leading to insufficient noise cancellation from elements with d = 7r₀ and beyond. Reducing this filter threshold to N_cut to 10⁻⁸ decreases the noise in the 6r₀ and 7r₀ elements of $U$ (see Fig. 12). Hence, if smaller elements need to be included in the generator when implementing STL-GME, one must be mindful of the numerical precision of the underlying dynamical solver.

Here, we summarize our modified protocol to build the STL-GME for a homogeneous 2D lattice. We note, however, that the same protocol can be used for an ND lattice with N ≥ 2. To construct the input dynamical matrix, C(t), we perform a single reference HEOM simulation and then exploit translational symmetry. Once we construct C(t), we follow the following workflow to implement the STL-GME for a 2D lattice:

1. Compute the generator $U(t)$ from reference C(t) using Eq. (8).

2. Identify the generator lifetime $τU$ by comparing the ‖L‖₂ error between TL-GME results and reference dynamics.

3. Identify the spatial distance cutoff $dU$⁠:

   • Keep the generator up to its lifetime, $U[N,N,t≤τU]$⁠.

   • Reshape the first two lattice indices into x and y coordinate indices such that any lattice index a can be rewritten as

     $a=i×Nx+j,$

     (D1)

     where N = N_x × N_y. This transforms the generator into a 5-rank tensor: $U[N,N,t]→U[Nx,Ny;Nx,Ny,t]$⁠.

   • To spatially truncate the generator for a homogeneous lattice,

     1. select an excitation starting at any lattice point, e.g., the origin of the 2D lattice, $U[Nx,Ny;0,0,t]$⁠,

     2. compute the distance between this initial excitation at coordinates (0, 0) and any measurement site (i, j) as

        $di,j=i2+j2,∀i,j,$

        (D2)

     3. discard all the elements where $di,j>dU$ by making them numerically 0,

     4. apply Eq. (C2) to this spatially truncated $U[Nx,Ny;0,0,t]$ to conserve the population, and

     5. employ translational symmetry to populate full $U[Nx,Ny;Nx,Ny,t]$ from $U[Nx,Ny;0,0,t]$ by varying initial excitation lattice coordinates.

   • Collapse the coordinate indices into lattice indices using Eq. (D1) such that $U[Nx,Ny;Nx,Ny,t]→U[N,N,t]$⁠.

   • Use the modified generator to propagate the dynamics using Eq. (6) (before $t<τU$⁠) and Eq. (7) after $t<τU$⁠.

   • Calculate the ‖L‖₂ error metric by comparing the exact dynamics with GME-predicted dynamics for different potential distance cutoffs to find $dU$⁠.

4. If $dU<N/2$⁠, the reference calculation is sufficiently large for spatial extension of the dynamics. To generate the dynamics of a system with extended size, perform the following steps:

   • Augment the generator’s dimension. Expand generator obtained from step (3c) $U[Nx,Ny;0,0,t]$ to $U[Mx,My;0,0,t]$ for all time points, where M_x ≥ N_x, M_y ≥ N_y and M = M_x × M_y, and assign 0 to all the new entries.

   • Similar to step (3d), use translational symmetry to populate $U[Mx,My;Mx,My,t]$ for all initial conditions from $U[Mx,My;0,0,t]$ and collapse the coordinate indices to lattice indices, yielding the expanded generator, $U[M,M,t]$⁠, for the M × M lattice system.

   • Similar to step 6 of  Appendix C, propagate the GME dynamics using $U[M,M,t]$⁠, which predicts the population correlation matrix C(t) for the M × M lattice.

For Figs. 8(b)–8(d), we simulate a 8 × 8 lattice to obtain the reference C(t) dynamics for the STL-GME extension. We find that the generator lifetime $τU=800$ fs and appropriate distance cutoff $dU=3r0$⁠, suggesting that the 8 × 8 lattice is sufficient as the reference simulation. We employ this 8 × 8 dynamics to predict the dynamics of a 10 × 10 lattice and confirm that our STL-GME dynamics agree with a separate exact HEOM simulation. In the main text [see Figs. 8(b)–8(d)], we employ our STL-GME protocol to augment the size to a 30 × 30 lattice and simulate over the first 100 ps.

Here, we discuss the spatial locality of the memory kernel in the time-local GME and how to truncate small elements that become negligible with increasing inter-site distance. For lattice problems with short-range interactions and local couplings (such as the dispersive Holstein lattice), space-locality is also present in the time-nonlocal description of memory. We demonstrate space-locality in the Transfer Tensor Method (TTM), which is the discrete analog of the time-nonlocal GME.¹¹¹ 

To demonstrate the performance of the SL-TNL-GME, we focus on the parameter regime of Fig. 1. Figure 14(a) confirms that the SL-TNL-GME can recover the reference simulation. Figure 14(b) shows that the characteristic memory distance in the time-nonlocal framework is $dK=5r0$⁠, which is larger than the analogous characteristic memory distance in the time-local framework, $dU=3r0$⁠. Similarly, we find the memory kernel lifetime in the time-nonlocal formulation is $τK=850$ fs, which is greater than the lifetime we identify in the time-local case, $τU=820$ fs. The latter result is consistent with previous findings in the context of biomolecular dynamics.⁶⁹ Hence, this example shows one can implement the idea of spatial memory truncation in the time-nonlocal formulation of reduced dynamics.