Lindblad dynamics from spatio-temporal correlation functions in nonintegrable spin-1 / 2 chains with different boundary conditions

,


I. INTRODUCTION
Quantum many-body systems out of equilibrium are a central topic of modern physics, and they have attracted increasing attention over recent years, both experimentally and theoretically [1][2][3][4][5].Key questions in this context are the emergence and properties of steady states in the limit of long times, but also the actual route to such states in the course of time [2][3][4][5].The understanding of these questions is of importance to isolated systems without any coupling to an environment and open systems with a weak or strong coupling to a bath, and it has witnessed a rather remarkable progress due to fresh concepts like eigenstate thermalization [6][7][8] and the typicality of random pure states [9][10][11][12][13][14][15][16][17] as well as due to the development of sophisticated numerical techniques [18,19].
In systems with a globally conserved quantity, a quite natural nonequilibrium process is given by transport [20].It is particularly a prime example of relevance to closed and open systems alike, in addition to the relevance of steady states and the relaxation to them.In isolated systems, a widely used approach is linear-response theory, which yields the well-known Kubo formula in terms of current autocorrelation functions [21].This theory can be formulated for density-density correlation functions as well, either in momentum and frequency space, or in the space and time domain.While linear response provides a clear-cut strategy, the concrete evaluation of correlation functions for specific models has turned out to be an analytical and numerical challenge, even for seemingly simple models of Heisenberg or Hubbard type in one dimension [20], with the most recent progress by generalized hydrodynamics [22,23].
In an open-system scenario, transport can be induced by the coupling of the system to baths at different temperatures or chemical potentials, which then usually yields a nonequilibrium steady state (NESS) with a constant current and characteristic density profile in the long-time limit [24][25][26][27].Such a scenario is often modelled by an equation of Lindblad form [28].While the derivation of this equation from a microscopic systembath model is a nontrivial task in practice [26,29], it is the most general version of a time-local master equation, which maps any density matrix to a density matrix again.In particular, it allows for an efficient numerical treatment by matrix-product states for quite large system sizes [24,[30][31][32], as dissipation reduces the unavoidable growth of entanglement as a function of time.
The two complementary approaches of closed and open systems have been the basis for a reliable identification of different dynamical phases, including the case of normal diffusion [20], but also for other transport [33][34][35][36][37][38][39][40][41][42][43][44] types ranging from ballistic to anomalous dynamics.Moreover, quantitative agreement of diffusion constants has been found in some cases [45][46][47][48].However, a clear correspondence between the dynamics in closed and open arXiv:2402.18177v2[cond-mat.stat-mech]7 Jun 2024 systems is still lacking [49][50][51].Thus, the main question of the present paper is: Can one predict the open-system dynamics based on the knowledge of the closed-system time evolution?While a general answer to this question appears to be hopeless, it might be possible for specific models and parameters at least.
To investigate this question, we rely here on a previous work [52], where spatio-temporal correlation functions have been suggested as a convenient ingredient in the case of weak driving and small system-bath coupling.Because the aforementioned work has focused on integrable systems and periodic boundary conditions, we intend to extend the analysis in three different directions: We (i) consider nonintegrable systems, (ii) take into account open boundary conditions and other bath-coupling geometries, see Fig. 1, and (iii) provide a comparison to time-evolving block decimation (TEBD).While we find that nonintegrability plays a minor role, the choice of the specific boundary conditions can be crucial, due to potentially nondecaying edge modes.Our large-scale numerical simulations suggest that a description based on closed-system correlation functions constitutes an useful alternative to existing state-of-the-art approaches.
Our paper is organized as follows.To begin with, we introduce in Sec.II the closed-system models studied and the spatio-temporal correlation functions.Then, in Sec.III, we discuss the concept of dynamical quantum typicality, and describe its implications for numerical and analytical purposes.Afterwards, in Sec.IV, we continue with the open-system setup, before we review the technique of stochastic unraveling in Sec.V.The subsequent Sec.VI is devoted to the central prediction used later, and its underlying assumptions.Eventually, in the next Sec.VII, we present our results.We conclude in Sec.VIII, and give additional information in the appendix.

II. CLOSED MODELS AND SPATIO-TEMPORAL CORRELATION FUNCTIONS
In this paper, we consider two different paradigmatic examples for quantum many-body models, which have attracted significant attention in the literature on, e.g., transport before.These two examples are nonintegrable modifications of the integrable spin-1/2 XXZ model in one dimension.The Hamiltonian of this model is given by [20] where S i r (i = x, y, z) are spin-1/2 operators at site r, N is the total number sites, J > 0 is the antiferromagnetic exchange coupling constant, and ∆ is the anisotropy in z direction.While the Hamiltonian in Eq. ( 1) is denoted for open boundary conditions, we will also use periodic boundary conditions, where the numbers of sites N is chosen to be even.The specific choice of boundary conditions will play an important role for the open setup, which will be discussed later in detail.
The spin-1/2 XXZ chain is well-known to be integrable for any value of ∆, and it has been in the focus of our previous work [52].In this work, we go beyond and include integrability-breaking perturbations, which then yield a more generic situation.As a first type of perturbation, we choose further interactions between next-to-nearest sites, which lead for open boundary conditions to and, although not shown explicitly, to a corresponding form for periodic boundary conditions.Here, ∆ ′ is the strength of the perturbation, and we will focus on the particular value ∆ ′ = 0.5, for which integrability is well broken.As a second type of perturbation, we choose a Zeeman term with a staggered magnetic field, where we will set the strength of the perturbation to the specific value B/J = 0.5, for the same reason as the one above.
Since the total magnetization S z = r S z r is strictly conserved for the models in Eqs.(3) and ( 4), [H, S z ] = 0, transport of spins is a meaningful question.Within the different approaches to transport in general, linear response theory is one of the main concepts.While this theory leads to the Kubo formula and current autocorrelation functions [21], it is also the basis for spatio-temporal correlation functions, Here, β = 1/T is the inverse temperature (measured in the units of the Boltzmann constant), and from now on we will consider the high-temperature limit β → 0, which still features nontrivial transport properties.The correlation functions in Eq. ( 5) measure the overlap of a timeevolved S z r (t) at some site r with an initial S z r ′ (0) at another site r ′ .In Fig. 2 (a), we illustrate the space-time dependence, as obtained numerically for the Hamiltonian H ′ in Eq. ( 3) with ∆ = 1.5, ∆ ′ = 0.5, N = 20, and periodic boundary conditions.Initially, there is a δ-function peak at the site r = r ′ and an uniform equilibrium background at other sites r ̸ = r ′ .The subsequent broadening of the δ-function peak corresponds to transport.
A convenient way to analyze the type of transport is provided by the spatial variance [53] Σ where the distribution C rr ′ (t) = 4⟨S z r (t)S z r ′ ⟩ eq is a normalized version of Eq. ( 5), r C rr ′ (t) = 1.For diffusive transport, Σ(t) ∝ t 1/2 .In Fig. 2 (b), we depict Σ(t) for the example in Fig. 2 (a).While it is clear that a ballistic growth Σ(t) ∝ t has to take place at short times, there is a diffusive growth Σ(t) ∝ t 1/2 at intermediate times, as expected due to the nonintegrability of the model [54].It is worth mentioning that the integrable model in Eq. (1) features a richer phase diagram, including superdiffusive behavior for ∆ = 1 and ballistic behavior for ∆ < 1.The different types of transport become also manifest in the quantity [53] which becomes constant in case of diffusion.Let us already mention here that the spatio-temporal correlation function ⟨S z r (t)S z r ′ (0)⟩ eq is not only a strategy to study transport in closed systems, but it may also be used to predict the build-up of a nonequilibrium steady state in open systems, where a bath is coupled at each edge.Investigating the quality of such a prediction is a central point of our paper.However, before we discuss this point in detail, we need to introduce further concepts in the following.

III. DYNAMICAL QUANTUM TYPICALITY
Next, let us discuss dynamical quantum typicality as one of the central concepts applied in this paper.On the one hand, this concept provides the basis for a numerical calculation of the spatio-temporal correlation functions ⟨S z r (t)S z r ′ (0)⟩ eq in closed systems of comparatively large  size.On the other hand, it constitutes a main ingredient to connect these correlation functions to the dynamics in open systems, as we will see in a later section.
Loosely speaking, the basic idea of typicality is that a single pure state can imitate the full statistical ensemble, on the level of the corresponding expectation values [9][10][11][12][13].To be precise, we introduce a pure state drawn at random (according to the Haar measure) from a Hilbert space of high dimension D, where {|j⟩} is an arbitrary orthonormal basis, and the real and imaginary part of the coefficients c j = a j + ib j result from a Gaussian probability distribution with zero mean and unit variance.For such a pure state, one then obtains the approximation [16,17] where the statistical error on the r.h.s. is exponentially small in system size, due to D = 2 N .This approximation can be rewritten as using the two auxiliary pure states |ψ(t)⟩ = e −iHt |ψ⟩ and |ϕ(t)⟩ = e −iHt S z r ′ |ψ⟩.The expression (10), just like analogous expressions for other observables, has turned out to be particularly useful for numerical simulations, since its evaluation requires the forward propagation of pure states in time.These propagations can be carried out efficiently in huge Hilbert spaces, which are orders of magnitude larger than the ones accessible by standard exact diagonalization [16,17].Note that the numerical data in Fig. 2 are also obtained in this way, yet on the basis of a slightly different and simpler expression, as explained in the following.
The simplification employs the fact that S z r , as well as n r = S z r + 1/2, are operators with the specific properties tr[S z r ] = 0 and n 2 r = n r .Using these two properties, and introducing the pure state then leads to expression [55] which involves a single pure state only.This expression is an obvious numerical benefit, but also provides a central analytical relation for later purposes in the context of open systems.

IV. OPEN SETUP AND LINDBLAD EQUATION
Now, let us turn to an open-system scenario, where we couple the system to an environment.We describe this scenario by the Lindblad equation, as the most general form of a time-local quantum master equation, which maps any density matrix to a density matrix again [28].The Lindblad equation ( 13) consists of a coherent part for the unitary time evolution w.r.t.
H and an incoherent damping term.This damping term is given by with non-negative rates α j , Lindblad operators L j , and the anticommutator {•, •}.Despite the generality of the Lindblad equation, its derivation is a challenging task for a specific microscopic model [26,29].
Here, we couple our system to two baths and choose the respective Lindblad operators [20] where γ is the system-bath coupling and µ is the driving strength.L 1 and L 2 are local operators, which act on a site B 1 and flip a spin up and down, respectively.L 3 and L 4 are corresponding operators at another site B 2 .
Our different choices of the bath-contact sites B 1 and B 2 are illustrated in Fig. 1.In case of periodic boundary conditions, we set B 1 = 1 and B 2 = N/2+1 in a distance N/2, see Fig. 1 (a).In case of open boundary conditions, we set B 1 = 1 at the left edge and B 2 = N at the right edge or, as an alternative, B 1 and B 2 close to the edges, see Figs. 1 (b) and (c).This alternative will be discussed in detail later.For all choices of B 1 and B 2 , the first (second) bath induces a net polarisation of order µ (−µ), leading to a steady state in the long-time limit with a characteristic density profile and a constant current.
In this open scenario, we are interested in the dynamics of local magnetization, as occurring at finite times and in the long-time limit.Thus, we study expectation values which depend on the parameters of the Hamiltonian H, but also on the two bath parameters µ and γ.As initial condition, we choose ρ(0) ∝ 1, which corresponds to the high-temperature limit β → 0 and a homogeneous profile of magnetization.

V. STOCHASTIC UNRAVELING OF THE LINDBLAD EQUATION
We aim at finding a solution of the Lindblad equation or an accurate approximation of the same.To this end, we rely on the concept of stochastic unraveling, which uses pure states |ψ⟩ rather than density matrices ρ [56,57].This concept consists of an alternating sequence of stochastic jumps with one of the Lindblad operators and, between the stochastic jumps, a deterministic time evolution w.r.t. an effective Hamiltonian, For our choice of Lindblad operators in Eqs. ( 15) - (18), this effective Hamiltonian takes on the form with the occupation number n r = S + r S − r = S z r + 1/2.In this work, we focus on the weak-driving case.For µ ≪ 1, H eff can be approximated as Hence, the time evolution w.r.t. to H eff becomes Therefore, the dynamics of a pure state is only generated by the closed system H, apart form the scalar damping term.This fact will be one of the main ingredients to connect the closed system and the weakly driven open system.For larger µ, the dynamics is more complicated and also involves the operators n B1 and n B2 .Since H eff is non-Hermitian, the norm of the pure state ψ(t) is not conserved and decays in the course of time, cf.
Eq. (23).Therefore, at some time t = τ , the condition ||ψ(t)⟩ || 2 > ε is first violated for a given ε, which is here drawn at random from a box distribution ]0, 1].At this point in time, a stochastic jump with one of the Lindblad operators takes place.The new and normalized pure state reads [57] |ψ where the specific jump is chosen with probability After the jump, the procedure continues with the next deterministic time evolution w.r.t.H eff .This sequence of stochastic jumps and deterministic evolutions leads to a particular trajectory |ψ T (t)⟩.The time-dependent density matrix according to the Lindblad equation can be eventually approximated by the average over different trajectories T .Thus, the expectation value reads where T max is a large enough number of trajectories.For T max → ∞, the approximation becomes an equality, and the stochastic unraveling can be used as an in principle exact numerical technique.Moreover, it will provide the basis for an analytical connection between closed-system and open-system dynamics.As already mentioned above, we consider ρ(0) ∝ 1 as the initial condition of the Lindblad equation, which is realized in the stochastic unraveling via a random pure state |ψ(0)⟩ of the form in Eq. (8).

VI. CONNECTION BETWEEN CLOSED AND OPEN SYSTEMS
A. Time evolution Now, we are ready to discuss our central prediction, which particularly connects the dynamics in closed and open systems.While this prediction has been presented in our previous work [52,58] in detail, we repeat here parts of the derivation, which will be of help to understand the results presented afterwards.
The derivation is based on the stochastic unraveling of the Lindblad equation.It is important to recall that we focus on the weak-driving case µ ≪ 1.In this case, the deterministic evolution w.r.t.H eff is unitary apart from a scalar damping term, cf.Eq. ( 23).As a consequence, when calculating the quantity this scalar cancels out.Because the initial condition is a random state of the form in Eq. ( 8), the first deterministic evolution is quite simple.Under unitary time evolution, a random states remains a random state with a uniform density profile, d r (t) = 0. Therefore, the first nontrivial event is the subsequent jump.
Without loss of generality, let us consider at some time t = τ a specific jump with one of the Lindblad operators, e.g., L 1 .Then, the resulting pure state reads which has exactly the same structure as the pure state in Eq. (11).Hence, we can employ dynamical quantum typicality and get where Θ(t) is the Heavyside function.By the use of the same arguments, we can obtain analogous relations for the remaining Lindblad operators L j , which then either involve ⟨S z r (t)S z B1 (0)⟩ eq or ⟨S z r (t)S z B2 (0)⟩ eq .Afterwards, averaging over all four jump possibilities yields (30) with jump probabilities p j = α j /4γ for a random state, cf.Eq. (25).Finally, by inserting the prefactors α j from Eqs. ( 15) -( 18), we end up with for the expectation value after the first jump.To proceed, a natural idea is to reuse the same line of reasoning for the second and all subsequent jumps.After the first jump, however, the pure state is different, since it has an inhomogeneous density profile with magnetization concentrated at the site of the bath contact.Thus, one has to wait until the injected magnetization has spread over a piece of the system.Clearly, such a waiting time requires a small enough value of the system-bath coupling γ, as illustrated in Fig. 3.This kind of equilibration is a central ingredient, and its impact will be scrutinized for specific models later.Assuming equilibration, we can iterate the arguments and obtain a superposition of the form where the amplitudes A j can be calculated from the jump probabilities in Eq. ( 25) and result as [52] A j = a j − dB1 (τ j − 0 + ) µ with Because the expression in Eq. ( 32) applies to a single sequence of jump times, (τ 1 , τ 2 , . ..), the final prediction is obtained by the average over trajectories with different jump times.Due to the scalar damping term in Eq. ( 23), these jump times are given by where ε j+1 is drawn at random from a box distribution ]0, 1] again.
In principle, the prediction in Eqs. ( 32) and ( 35) can be calculated analytically for a specific model.However, the closed-system correlation functions ⟨S z r (t)S z B1 (0)⟩ eq and ⟨S z r (t)S z B2 (0)⟩ eq are often available only numerically, such that the prediction has to be calculated numerically as well.

B. Injected magnetization
While Eqs. (32) and (35) allow to predict the dynamics of magnetization at finite times and in the limit of long times, a similar expression can derived for the respective currents in the steady state.To this end, let us consider the magnetization injected by the first bath, which can be predicted as [52]  with δ being just a notation for "injected" and which is slightly simpler than Eq.(32).Since in the steady state all currents are the same, it is sufficient to know ⟨j B1 ⟩, which can expressed as Here, f = 2 for periodic boundary conditions (flow to the right and left of the bath) and f = 1 for open boundary conditions (flow only to the right of the bath).
With the knowledge of the steady-state current, it is also possible to predict the diffusion constant via [20] for some site r in the bulk.Eq.TEBD Eq. TEBD Eq. TEBD FIG. 7. Similar comparison as the one in Fig. 6, but now, instead of stochastic unraveling, between the prediction (Eq.) and time-evolving block decimation (TEBD).

A. Next-to-nearest-neighbor interactions and periodic boundary conditions
Finally, we turn to our numerical simulations, where the central goal is to analyze the quality of the prediction in Eqs.(32) and (35) for various situations.To start with, we investigate the spin-1/2 XXZ chain with interactions between next-to-nearest neighbors, Eq. ( 3), and the case of periodic boundary conditions, Fig. 1 (a).Afterwards, we additionally study other perturbations and the case of open boundary conditions with different bath-coupling geometries.
Because a main ingredient of the prediction has been the equilibration of the injected magnetization, we first focus on this assumption.To this end, we numerically calculate in Fig. 3 (a) the equal-site correlation function ⟨S z r (t)S z r (0)⟩ eq for ∆ ′ = 0.5, N = 20, and an arbitrary site r due to periodic boundary conditions.Apparently, for all ∆ depicted, this function starts with the initial value 1/4, decays substantially on a time scale t R J ≈ 5, and then approaches the equilibration value 1/(4N ) in the limit of long times.By comparing t R to 1/(2γ) from Eq. ( 36), we identify γ/J = 0.1 as a reasonable choice for the system-bath coupling, which we fix from now on for a fair comparison, together with the choice µ = 0.1 to ensure weak driving.
For the value ∆ = 0.5, we depict in Fig. 4 (a) the time evolution of magnetization in the open system.Here, the prediction is carried out for ≈ 10, 000 different sequences of jump times, which already yields smooth curves.The full stochastic unraveling, without any assumption, has been evaluated on clusters and turns out to require as many as ≈ 200, 000 (or more) sequences for a comparable smoothness [59].Despite residual statistical fluctuations, the agreement is almost perfect and clearly visible from short until long times.This convincing agreement can also be seen for the steady-state profile in Fig. 4 (b).Compared to previous results for integrable systems [52], the degree of agreement turns out to be similar.This finding indicates that nonintegrability is not required for the prediction to hold, which is consistent with the fact that such an assumption does not enter the derivation.
In Fig. 4 (c), we also show the injected magnetization, which grows linearly with time.From the slope, and the slope of the steady-state density profile in Fig. 4 (b), we can extract a diffusion constant via Eq.( 41).The value D/J ≈ 2.9 agrees well with other values in the literature, such as D/J ≈ 3.1 in the closed system [60,61], and serves as a further sanity check of our approach.
To ensure that our results do not depend on the specific choice of parameters, we redo in Fig. 5 the calculation for other ∆ ̸ = 0.5.The overall agreement is apparently the same.

B. Open boundary conditions and different bath-coupling geometries
Now, we move forward to the case of open boundary conditions with a standard bath coupling at the edges, as sketched in Fig. 1 (b).For such a situation, it is also possible to obtain the solution of the Lindblad equation by TEBD [19,62], as a state-of-the-art technique in this context.For the same Hamiltonian and parameters as before, we depict in Fig. 6 the numerical result from TEBD and additionally compare to the exact stochastic- unraveling procedure.Even though statistical fluctuations are again visible, both approaches coincide for all times and ∆ depicted.This agreement particularly confirms the correctness of our numerics.
When we compare TEBD to the actual prediction in Fig. 7, the agreement is less convincing for long times and becomes worse for larger values of ∆.To understand the origin of the disagreement, we test the assumption of equilibration, by calculating in Fig. 3 (b) the equal-site correlation function ⟨S z r (t)S z r (0)⟩ eq at the left-edge site r = 1.In contrast to periodic boundary conditions, this function decays slower and develops long-time tails for large ∆, which can be traced back to nondecaying edge modes occurring for ∆ > 1 [63,64].As the occurrence of long-time tails in Fig. 3 (b) seems to correlate with the degree of disagreement in Fig. 7, we can identify the breakdown of the assumption as the origin.
For ∆ = 0.5, where the assumption is fulfilled best, we increase next the system size to N = 34 ≫ 20.For such a system size, the Hilbert-space dimension becomes huge and stochastic unraveling is not feasible any more, due to the many trajectories required.However, the prediction can still be carried out, since the two correlation functions ⟨S z r (t)S z B1 (0)⟩ eq and ⟨S z r (t)S z B2 (0)⟩ eq need to be calculated only once.In particular, this calculation is possible by the use of dynamical quantum typicality and supercomputers.In Fig. 8, we depict the corresponding prediction and compare to the TEBD solution, which up to times tJ ≈ 50 does not depend on the bond dimension used.(A convergence analysis of TEBD can be found in the appendix.)The convincing agreement sup- ports that our prediction is a useful alternative for large system sizes, which are usually only accessible by matrix product states.
Unfortunately, our assumption is not always satisfied for open boundary conditions, as discussed above.Hence, we explore possibilities to circumvent this problem.To this end, we consider the slightly different bath-coupling geometry in Fig. 1 (c), where the Lindblad operators are not located exactly at the edges, but close to them.As depicted in Fig. 9 (a), the equal-site correlation function ⟨S z r (t)S z r (0)⟩ eq tends to decay stronger when the site r is moved away from the left edge r = 1, indicating a higher degree of equilibration.And indeed, for the so far worst case ∆ = 1.5, the prediction for the open-system dynamics in Fig. 9 (b) agrees substantially better with the full stochastic unraveling.This observation supports the usefulness of other bath-coupling geometries, which have attracted less attention yet.

C. Staggered field
Eventually, we also study other perturbations and turn to the spin-1/2 XXZ chain with a staggered field, Eq. ( 4), where we focus on ∆ = 1.0,B/J = 0.5, and N = 20.In Fig. 10, we summarize our numerical results.Apparently, the situation is overall similar.The equal-site correlation function ⟨S z r (t)S z r (0)⟩ eq in Fig. 10 (a) behaves differently for periodic and open boundary conditions.Consistently, the prediction for the open-system dynamics agrees well with numerics for the periodic-boundaries case in Fig. 10 (b), while deviations are visible for the open-boundaries case in Fig. 10 (c).We have checked that the situation remains the same for other parameters of ∆, even though not explicitly shown here.

VIII. CONCLUSION
To summarize, we have studied the Lindblad equation as a central approach to boundary-driven magnetization transport in spin-1/2 chains.Our main motivation has been to understand to what extent the dynamics in the open system, at finite times and in the limit of long times, can be predicted on the basis of the dynamics in the closed system.To this end, we have followed the idea of a previous work [52], which has suggested a prediction in terms of spatio-temporal correlation functions, Eqs. ( 32) and (35), given the case of weak driving and small system-bath coupling.While this work was focused on integrable systems and periodic boundary conditions, we have substantially extended the analysis in the current work by going in three different directions: We have (i) considered nonintegrable systems, (ii) included open boundary conditions and other bath-coupling geometries, and (iii) provided a comparison with time-evolving block decimation.
We have found that nonintegrability plays a minor role, since the quality of the prediction is comparable to the one for integrable systems.This observation is consistent with the fact that nonintegrability does not enter as an assumption in the derivation.In contrast, the choice of the specific boundary conditions has turned out to be of relevance.For periodic boundary conditions, on the one hand, prediction and numerical simulations have agreed convincingly, for all models and parameters considered here.For open boundary conditions, on the other hand, we have observed some disagreement in particular cases, which we have traced back to slowly decaying edge modes and thus a breakdown of the equilibration assumption underlying the prediction.In this context, it is important to note that the validity of the assumption can be checked in advance and does not require the comparison to other methods.To circumvent such edge modes, we have also explored other bath-coupling geometries, where the Lindblad operators are not acting exactly at the boundary sites, but still close to them.
For parameters, where the assumption is well fulfilled also for open boundary conditions, we have demonstrated that the prediction yields an accurate description and can be carried out for comparatively large system sizes, which are usually accessible by matrix-product-states methods only.From a less practical but more physical perspective, we have thus shown a kind of one-to-one correspondence between the time evolution in open and closed systems, at least for the models considered by us.
Promising future directions of future research include quasi-1D lattices, finite temperatures, energy transport, fermionic models, and disorder.Another interesting avenue is to explore in which cases the slowly decaying edge modes of the closed system can be enhanced in the open system [65].no.101060162, and the Packard Foundation through a Packard Fellowship in Science and Engineering.We gratefully acknowledge the Gauss Centre for Supercomputing e. V. for funding this project by providing computing time on the GCS Supercomputer JUWELS6 at Jülich Supercomputing Centre (JSC).S. N. was supported by QuantERA grants QuSiED and T-NiSQ, by MVZI, QuantERA II JTC 2021.TEBD calculations were performed on the supercomputer Vega at the Institute of Infor-mation Science (IZUM) in Maribor, Slovenia.

Appendix: Convergence of the TEBD method
In the main text, we have shown in Fig. 8 numerical data from TEBD and stated that, up to times tJ ≈ 50, it does not depend on the bond dimension χ used.To further substantiate this statement, we depict in Fig. 11 the same data for different χ and various sites, close to the edges and in the bulk.While the data is particularly well converged close to the edges, some oscillations can be seen in the bulk, where the time evolution is kind of close to unitary.

Appendix: Correlation functions for open boundary conditions
In Fig. 3 (b), we have shown the equal-site correlation function ⟨S z r (t)S z r (0)⟩ eq for the open-boundaries case and different ∆, where we have focused on a single system size N = 20.To demonstrate that the temporal decay does not depend significantly on system size, we additionally depict in Fig. 12 (a) numerical data for N = 24.
For completeness, Fig. 12 (b) shows the full time-space dependence of the correlation functions ⟨S z r (t)S z r ′ (0)⟩ eq for N = 34 and ∆ = 0.5, which have been used for the prediction in Fig. 8.

FIG. 1 .
FIG. 1. Sketch of the three different geometries of the bath coupling.(a) Periodic boundary conditions.(b) and (c) Open boundary conditions with a bath coupling located exactly at the edges or close to the edges of the system.Note that open boundary conditions should not be confused with open-system scenario.

FIG. 4 .
FIG. 4. Open-system dynamics for the model H ′ in Eq. (3), obtained numerically for ∆ = 0.5, ∆ ′ = 0.5, N = 20, periodic boundary conditions, small coupling γ/J = 0.1, and weak driving µ = 0.1.Exact results from the full stochastic unraveling (data) are compared to the prediction (Eq.), which is based on spatio-temporal correlation functions in the closed system.(a) Time evolution of the local magnetization ⟨S z r (t)⟩ for different sites r.(b) Site dependence of the steady state at tJ = 50.(c) Magnetization injected by the first bath as a function of time.

FIG. 8 .
FIG.8.Similar comparison as the one in Fig.7(a), but now for system size N = 34 ≫ 20, where stochastic unraveling is not possible any more.

FIG. 9 .
FIG. 9. On a bath coupling close to the edges.(a) Temporal decay of the autocorrelation function ⟨S z r (t)S z r (0)⟩eq for different sites r, as obtained numerically for the model H ′ in Eq. (3) with ∆ = 1.5, ∆ ′ = 0.5, N = 20, and open boundary conditions.(b) Open-system dynamics for a bath coupling at sites B1 = 3 and B2 = 17.In comparison to the data in Fig.7 (c), the agreement between the prediction (Eq.) and stochastic unraveling (data) is better.

FIG. 11 .
FIG. 11.Convergence analysis.The data in Fig. 8 is depicted again for different sites: (a) r = 17, (b) r = 13, (c) r = 9, and (d) r = 5.But now, the prediction (Eq.) is compared to data from time-evolving block decimation (TEBD) for various bond dimensions χ.While r = 5 is close to the edge, r = 17 lies in the bulk.

FIG. 12 .
FIG. 12. (a) Dynamics of the autocorrelation function ⟨S z 1 (t)S z 1 (0)⟩eq for open boundary conditions, as depicted in Fig. 3 (b) for different ∆, but now for two system sizes N = 20 and N = 24.(b) Time-space density plot of the correlation functions ⟨S z r (t)S z 1 (0)⟩eq for N = 34 and ∆ = 0.5, which are used for the prediction in Fig. 8.