Continuous and time-discrete non-Markovian system-reservoir interactions: Dissipative coherent quantum feedback in Liouville space

Based on tensor network realizations of path integrals reducing exponential memory scaling to polynomial efficiency and a Liouville space implementation of a time-discrete quantum memory, we investigate a quantum system simultaneously exposed to two structured reservoirs. For this purpose, we employ a numerically exact quasi-2D tensor network combining both diagonal and off-diagonal system-reservoir interactions with a twofold memory for continuous and discrete retardation effects. As a possible example, we study the non-Markovian dynamical interplay between discrete photonic feedback and structured acoustic phonon modes, resulting in emerging inter-reservoir correlations and long-living population trapping within an initially-excited two-level system.

In this work, we present a numerically exact tensor network-based approach allowing for the first time to describe two non-Markovian processes simultaneously, namely continuous and discrete retardation effects where interactions with both diagonal and off-diagonal system couplings are taken into account, i.e., couplings without or with energy exchange between system and reservoir. Recently established matrix product state (MPS) techniques to implement a time-discrete memory were aimed at quantum systems where decoherence and dephasing effects are not a key factor, and thus were based on a wave function ansatz to solve the quantum stochastic Schrödinger equation [27][28][29]. We extend this approach to a time bin-based density matrix description in Liouville space to include Markovian and non-Markovian decoherence effects. In a second step, it is combined with a tensor network-based real-time path integral method to describe interactions with a continuous structured reservoir [30][31][32][33][34], resulting in a quasi-2D tensor network formalism [35][36][37][38]. This architecture enables non-Markovian simulations of quantum systems coupled to two structured reservoirs, see Fig. 1, maintaining the relevant entanglement information and capturing both diagonal and off-diagonal system-reservoir interactions on equal footing. Possible applications include setups of waveguide-QED with dephasing [25,[39][40][41], e.g. realized by additional decay channels, or multiple spatially separated solid-state quantum emitters initially prepared in a dark state and interacting with their environment [42]. Here, we specifically consider a twolevel quantum system coupled to a structured reservoir of independent oscillators and subject to time-discrete coherent quantum feedback, extending the paradigm of the Spin-Boson model to the feedback realm [43]. We demonstrate that non-Markovian interplay between relaxation and decoherence processes results in a dynamical protection against destructive interference and thereby enables population trapping. This expands upon the widely-discussed localized phase stabilization in the spinboson model [31,43] from an incoherent feedback-induced perspective, replacing the coherent driving with another structured reservoir.
The paper is organized as follows: In Sec. I, a recently established MPS-based path integral implementation for continuous reservoirs is discussed. Afterwards, we introduce an MPS implementation of a time-discrete quantum memory in Liouville space in Sec. II. In Sec. III, we combine the two algorithms to form a quasi-2D tensor network, enabling numerically exact calculations of two non-Markovian system-reservoir interactions, before demonstrating its capabilities in Sec. IV, where we find a dynamical protection of coherence in the presence of two non-Markovian reservoirs. Lastly, we summarize our findings in Sec. V.

I. PATH INTEGRAL FORMULATION FOR CONTINUOUS RESERVOIRS
We start with the implementation and evolution of a system subjected to a continuous harmonic reservoir. For a numerically exact description, our theoretical approach is based on a real-time path integral formulation [9][10][11][12][13][14][15][16][17]. In recent breakthroughs, path integrals have been implemented in a tensor network approach based on MPS [31,32], allowing to solve non-Markovian dynamics and providing an efficient representation using highdimensional tensors with restricted correlations [31]. In the following, we briefly recapitulate the algorithm introduced by Strathearn et al. [30,31,34], which is employed as a part of our solution to multiple non-Markovian system-reservoir interactions.
Our goal is to employ path integrals for a numerically exact solution of the von-Neumann equation given a Hamiltonian H(t) describing a time-dependent systemreservoir interaction [44,45], with ρ(t) denoting the density matrix and L(t) the Liouvillian superoperator. As a main challenge, the evaluation of path integrals becomes increasingly expensive over time, since the history of all preceding paths at times 0, . . . , t n−1 must be taken into account for the calculation of the current time step t n . However, in case the system-reservoir correlations are finite in time, the augmented density tensor scheme can be introduced for improved numerical accessibility [16,17]. Exploiting the finite reservoir memory length, only the last n c time steps are taken into account for the calculation of the current path. This treatment is known as the finite memory approximation and results in the augmented density tensor representation as a solution to the system part of Eq. (1) with traced out reservoir contributions, which reads at and constitutes a discrete path integral formulation where indices i (′) n contain the left (right) configuration of the system at time t n = n∆t. The field transformation matrix M inin−1 e.g. accounts for the action of an external driving field Ω 0 . For the case Ω 0 = 0 considered below, it takes the simple form M inin−1 = 1δ in,in−1 . The influence functional is given by with and φ(τ − τ ′ ) the reservoir autocorrelation function [10,11,43]. Using an improved finite memory approximation, for n − m ≡ n c all former paths up to t nc = n c ∆t are additionally incorporated in the integration, i.e. η nc := η n−m + n−nc−1 k=1 η n−k [30]. Under this approximation, it is possible to restate the augmented density tensor and its time evolution efficiently as a tensor network [31]: First, Eq. (2) is mapped to a vector ρ jn in Liouville space, with I(j n , j m ) := jn−1M jnjn−1 exp (S jnjm ). Here, left and right system indices i k , i ′ k have been combined to a single index j k for each time step, resulting in Liouville space representationsM jnjn−1 andS jnjm of the field transformation matrix and the influence functional, respectively. Afterwards, the augmented density tensor is rewritten as an MPS, storing the present and up to n c − 1 past states in individual tensors with the oldest state located at the left end of the MPS. In this representation, tensor compression by consecutive applications of the singular value decomposition [46] reduces the memory requirements to polynomial rather than exponential scaling with respect to n c [31].
The time evolution is carried out by a network of matrix product operators (MPOs), shown schematically in Fig. 2(a) (dark grey shapes). During the first time step, the initial system state ρ j0 (0) (red shape) is contracted with the first MPO in the network [dashed frame in Fig. 2(a)]. As a result, the system state is updated and the preceding path is stored to its left, increasing the length of the MPS by one. Fig. 2(b) shows the MPS after completion of the first time step. Once step n = n c is reached, the oldest path in the MPS is summed over by application of a delta tensor [semicircular shape in Fig. 2(a)], corresponding to the improved finite memory approximation [30,31]. At this stage, the MPS length is fixed for the rest of the time evolution. Moreover, for time-independent problems, apart from the index nomenclature the structure of the MPO remains unchanged for all time steps n ≥ n c , resulting in an additional performance gain.
To provide an example of a continuous reservoir of noninteracting harmonic oscillators, we consider the Hamiltonian [47] corresponding to diagonal system coupling without inducing system transitions, with system operators σ ij = |i j|, bosonic annihilation (creation) operators b ( †) q of reservoir modes with frequency ω q = c s |q|, c s the sound velocity, and a mode q-dependent system-reservoir coupling amplitude g q . The corresponding correlation function reads with temperature T and k B the Boltzmann constant.
In the following, we choose a generic coupling element e.g. describing acoustic bulk phonons interacting with a quantum emitter [48]. As a benchmark of the tensor network-based path integral implementation, we first calculate the analytically solvable Independent Boson model, consisting of a single two-level emitter subjected to pure dephasing by a structured harmonic reservoir, as described by Eq. (6). In Fig. 2(c), we prepare the initial polarization at ρ 01 (0) = 0.5i and calculate the resulting dynamics at varying temperatures T (solid lines). The corresponding analytical solution (dashed grey lines) is given by [49] ρ 01 (t) = exp (8) exhibiting excellent agreement with the numerical results at all considered temperatures. In addition, the employed method features very high performance, enabling reservoir memory depths of n c = 100 and beyond. To exemplify the capabilities of the tensor network implementation, we calculate the time evolution dynamics of the Spin-Boson model, corresponding to Eq. (6) with an additional continuous driving field term at amplitude Ω 0 , Fig. 2(d) shows the resulting dynamics at parameters n c = 100, Ω 0 = 0.5 ps −1 and T = 77 K for the memory, driving field and temperature, respectively.

II. TIME-DISCRETE MEMORY IN LIOUVILLE SPACE
As a second non-Markovian reservoir, we consider a discrete time-bin based quantum memory. Recently established implementations rely on an MPS-based wave function ansatz to compute the quantum stochastic Schrödinger equation [27][28][29]. However, this formulation is inherently incompatible with the previously introduced path integral formulation. As a solution to this problem, we present an MPS implementation of a time-discrete quantum memory in Liouville space. Here, the dynamics of the system density matrix is prescribed by a Liouvillian superoperator [see Eq. (1)], with a Hamiltonian H D containing the time-delayed system-reservoir coupling, such that interactions occurring at time t couple back into the system and affect its state at a subsequent time t + τ , with τ the retardation time. Such a time-discrete coupling e.g. arises in a two-level emitter with states |0 , |1 at an energy difference ω 0 , placed in front of a mirror with a round trip time τ . The corresponding Hamiltonian reads describing off-diagonal system coupling leading to energy exchange between system and reservoir and system phase relaxation, with system operators σ ij = |i j|, bosonic annihilation (creation) operators r ( †) k of photon modes with frequency ω k = ck, c the speed of light, and a constant electron-photon coupling amplitude Γ. For an efficient evaluation, the dynamics imposed by L(t) is translated in a time bin-based MPS formalism [27][28][29]37] which maintains the relevant system-environment correlations, scaling with τ . In case of additional phenomenological dissipative channels, the Liouvillian can be extended by the standard Lindblad operator [44,45].
For the MPS implementation of the time-ordered Liouvillian, we start from the formal solution of the system part of Eq. (1) for the density matrix, with T the time-ordering operator. For an MPS-based approach and in analogy to the time-discrete path integral formulation, we restate Eq. (10) in a time-discrete basis, which reads at time t N = N ∆t with at time discretization ∆t and with time-bin normalized operators For the MPS evolution of the density matrix during each time step n, the discrete Liouvillian time step operator L(n, n − 1) is approximated as a tenth order series expansion, i.e., Figs. 3(a)-(e) show the tensor network scheme for the implementation of the time-discrete memory. The square red tensor in Fig. 3(a) contains the system density matrix at the initial time t = 0. To consider a time-discrete memory, here n d = 4 circular tensors to its left store the reservoir state in Liouville space at preceding times, with the oldest state located on the left end of the MPS (blue). The reservoir states for all future time steps are initialized to the right of the system bin, containing full reservoir entanglement e.g. at finite temperature. In the following we assume an initial vacuum state. Therefore, during each time step a new empty reservoir bin (orange) is added to the MPS from the right, representing the present reservoir state [see Fig. 3(a)]. The memory loop realization explained in detail below introduces the retardation time τ = n d ∆t by the number of initial memory bins n d .
The first step of the time evolution is carried out as follows: By applications of the singular value decomposition algorithm [46], the first memory bin (blue) is pushed to the left of the system (red) while maintaining relevant entanglement information in the swapping procedure [see Fig. 3(a)]. The Liouvillian operator L(1, 0) for the first time step is then applied to the system bin, current memory bin and present reservoir bin, as shown in Fig. 3(b). Afterwards, the processed memory bin (grey) is swapped back to its original position and stored for the rest of the time evolution. The updated present reservoir bin (green) is pushed to the left, taking the role of a memory bin [see Fig. 3(c),(d)]. Fig. 3(d) shows the MPS after completion of the first time step. The second time step is carried out in the same fashion, as shown in Fig. 3(e). After completion of n d time steps, all initial memory bins have been processed. At step n d + 1, the reservoir bin modified during the first time step [green bin in Fig. 3(d)] becomes the current memory bin, containing information of a previous system state and setting off reservoir-induced memory effects in the system in complete agreement with the time-ordered problem.
As a first benchmark for the presented time-discrete quantum memory in Liouville space, we calculate the system dynamics imposed by Eq. (9). Fig. 3(f) shows the unfolding emitter population dynamics at a feedback time τ = 3.0 ps and Γ = 31.6 ps −1 calculated using the MPS implementation (solid blue line) and compared to its analytical solution up to t = 3τ (dashed grey line). The latter is given by σ 11 (t) = | σ 10 (t) | 2 only valid in the single-excitation regime with (14) as calculated in [50][51][52] and shows excellent agreement with the numerical result. In this regime, i.e. Γτ ≫ 1, the delay in the amplitude governs the dynamics and leads to re-excitations at multiples of the round trip time τ . The phase of the amplitude φ = ω 0 τ , however, loses importance in the first τ -intervals due to a stronger decay of the mixing terms in the absolute square of Eq. (14). The advantage of our implementation of quantum feedback in Liouville space becomes evident if the impact of phase destroying processes is in question. Up until now, this impact has only been investigated for a special case, finding the emergence of an Ornstein-Uhlenbeck process during the first τ -intervals [52]. These results have been obtained via analytical calculations, limiting the investigation to a small number of feedback intervals. Steadystate scenarios, however, are out of reach in this case as the evaluation of the phase-noise kernels must be done analytically. In our method, these limits have been overcome. Due to the here presented Liouville architecture, additional Lindblad-based dissipation can be easily implemented without increased numerical expense. In Fig. 3(g), we present the dynamics of a decaying, intially excited two-level emitter under the influence of quantum coherent feedback and additional phenomenological dephasing at rate γ, realized by adding a Lindblad dissipator to the Liouvillian [Eq. (1)] and H = H D , [44,45] with a redefined system operator in full configuration space,σ 11 = 1 D σ 11 1 D , including the time-discrete reservoir basis 1 D = dk ∞ n=0 |{n k } {n k }|. The time trace shows long time calculations of the emitter population, comparing the cases γ = 0 (blue line) and γ = 0.5 ps −1 (orange line). Most importantly, we see that the pure dephasing process becomes important only after the feedback signal re-excitates the emitter and the stabilization of the incoming and outgoing phase comes into play. In the presence of an additional pure dephasing γ = 0, the initial decay process is unchanged but the re-excitation becomes less efficient until only incoherent re-excitation takes place, leading to a faster decay to zero without population trapping, regardless of the choice of φ.
This important result sheds light on the robustness of quantum feedback processes in the presence of additional Markovian dissipation channels. As expected, additional Markovian decoherence leads to a washing out of the signal since a loss of quantum feedback-induced coherence is inevitable. However, this does not have to be the case in the presence of an additional non-Markovian dissipation channel, which we discuss in the following.

III. QUASI-2D TENSOR NETWORK
As a next step, we expand upon the MPS architecture for time-discrete quantum memory in Liouville space by combining it with the previously discussed tensor network-based path integral implementation for continuous harmonic reservoirs, resulting in a quasi-2D tensor network. The technical connection of the networks via link indices which store arising entanglement information enables the numerically exact description of correlation buildup in between the reservoirs. When considering a scenario involving two non-Markovian reservoirs not isolated from each other, such inter-reservoir correlations may have fundamental impact on the system dynamics, therefore prohibiting a strict truncation of the arising inter-reservoir entanglement, e.g. in the form of a low Schmidt value cutoff precision d cut . As a result of not only two system-reservoir interactions but additionally arising reservoir-reservoir entanglement, the overall grade of entanglement in the system rises intensively with respect to the twofold single reservoir case. On the other hand, in setups where two non-Markovian reservoirs are present but do not crucially interact, e.g. via dynamical decoupling, a much more restrictive truncation is possible without cost of accuracy. For the presented results, we have employed a high Schmidt value cutoff precision d cut = 10 −12 , such that no relevant entanglement information is lost during the time evolution. The dynamical interplay between the two reservoirs with the system and with each other poses an immense numerical challenge and strongly limits the accessible memory depths in the here considered system: While the two presented algorithms by themselves enable simulations of a single reservoir with deep memories, their combination is accompanied by limitations due to the arising interreservoir entanglement. As a result, the combined number of memory bins in the quasi-2D network is limited to n c + n d < 20 for our model of choice, as is the case in traditional path integral implementations for a single reservoir [11][12][13][14][15]43]. However, we stress once more that this limitation is a natural consequence of the high grade of entanglement in between the two reservoirs and the system. With the presented quasi-2D network architecture, we take first steps to unravel the mostly unexplored field of multiple interacting non-Markovian reservoirs by explicitly considering memory-enabled information backflow in between them.
The construction of the quasi-2D network is sketched in Fig. 4: The two tensor networks for the continuous [Figs. 4(a),(b)] and time-discrete reservoirs [ Fig. 4(c)] are connected to each other via the common tensor representing the current system state in both MPS algorithms (red shape). The system state tensor is employed to act as a junction connecting the two reservoirs [dashed circles in

Figs. 4(b),(c)]
and thereby enables the buildup and storage of inter-reservoir correlations in the connecting link indices. The time evolution of the quasi-2D network is carried out as follows: During each time step, the system is first evolved under the influence of the continuous reservoir by a single contraction of the network, as shown in Fig. 4(a) (dashed frame). Afterwards, the new current system state [red shape in Fig. 4(b)] is subjected to the second tensor network algorithm accounting for the time-discrete reservoir [ Fig. 4(c)]. In consequence, the quasi-2D network stores the history of both interactions, maintaining crucial entanglement information and enabling the calculation of two dynamically interacting time-delayed processes. As an example application for the quasi-2D network, in the following we investigate the interplay of off-diagonal coherent quantum feedback and a diagonal reservoir of independent oscillators. As illustrated below, for certain memory depths and initial states, this results in dynamical protection of coherent quantum feedback properties in the open system.

IV. MEMORY-INDUCED DYNAMICAL POPULATION TRAPPING
Coherent quantum feedback mechanisms exhibit a rich variety of non-Markovian phenomena [51,[53][54][55][56][57][58], e.g. enabling coherent population trapping [28,42,[59][60][61][62], Ornstein-Uhlenbeck-type events in the presence of white noise [52], and formation of large entangled photon states [63]. However, so far these effects have not been explored in the presence of additional non-Markovian decoherence or dissipative channels. To investigate the impact of dephasing on feedback-induced decoherence, we consider a two-level emitter placed in front of a mirror with a round trip time τ , taking the role of a time-discrete reservoir (see Fig. 1) described by Eq. (9). The photoninduced feedback imprints a time-delayed coherence in the form of a feedback phase ϕ = ω 0 τ /(2π) on the system, critically influencing its dynamics. It is given by the delay time τ and the transition frequency ω 0 of the electronic coherence operator σ 12 . To study a pronounced quantum optical effect, we consider the case of an initially excited two-level system where coherent population trapping occurs as a result of a bound state in continuum at feedback phases ϕ ∈ Z [28,42,52,[59][60][61][62]. The unfolding dynamics are evaluated for parameters Γ = 0.9 ps −1 , τ = 1.2 ps and n d = 4, resulting in a time discretization ∆t = 0.3 ps, e.g. typical for semiconductor quantum dot based devices [58,61,62]. Moreover, we get Γτ ≈ 1.1, corresponding to the strong non-Markovian regime [27,52,64]. The employed series expansion of the Liouvillian up to tenth order [see Eq. (13)] justifies this coarse time discretization, making n d = 4 time bins sufficient for our investigation while resulting in convergent results (see Appendix A).
To illustrate the power of our method, we compare the cases of Markovian and non-Markovian dephasing introduced by an additional diagonal system-reservoir coupling, representing the continuous reservoir in Fig. 1. As a first step, we calculate the system dynamics given by Eq. (9) in the presence of phenomenological dephasing at rate γ, introduced by the Lindblad dissipator stated in Eq. (15). Fig. 5 shows resulting excited state population dynamics at varying γ and feedback phases ϕ. At γ = 0 there exists a periodic ideal feedback phase ϕ ∈ Z such that the system decouples from its environment by constructive interference, resulting in coherent population trapping (solid blue line in Fig. 5). Choosing a nonzero dephasing γ = 0.001 ps −1 has no impact on the population dynamics until feedback sets in, since the radiative decay is frequency-independent until t = τ (dashed blue line). Thereafter, phenomenological dephasing destroys the phase interference and with it the trapping mechanism, resulting in an asymptotic decline of occupation to zero. At a nonideal feedback phase, here ϕ = 1.17, and no dephasing, destructive interference leads to an asymptotic decline towards zero as well (solid orange line). The decay is further accelerated by setting γ > 0 (dashed orange line), since any phenomenological decoherence attacking the phase relation ϕ results in faster decay. In conclusion, choosing a feedback phase ϕ / ∈ Z without a structured phonon reservoir inevitably results in asymptotic population decay via thermalization, and a Lindblad formulation of decoherence never preserves quantum correlations between the reservoir and system states.
As a next step, we show this is not necessarily the case if the decoherence process itself is the result of a non-Markovian reservoir interaction. Using the quasi-2D tensor network, we calculate the emitter dynamics imposed by H D + H C [Eqs. (6), (9)] at n c = 4 and leaving all remaining parameters unchanged. Fig. 6 shows resulting population dynamics at different temperatures. For the chosen parameters, we find population trapping for the non-ideal feedback phase ϕ = 1.17, i.e. ϕ / ∈ Z, at T = 4 K (solid line). Time-delayed excitation backflow from the continuous reservoir of oscillators to the system enables a decoupling from destructive interference with the time-discrete photon environment. This information backflow results in correlation buildup and information exchange between the reservoirs, dynamically protecting feedback-induced coherence in the system for long times. As long as diffusion processes at finite temperature take place on a comparable time scale as the coherent feedback dynamics, we always find dynamical population trapping by tuning of ϕ after a typical excitation backflow time. At higher temperatures, it is to be expected that thermal properties of the phonon reservoir start to dominate the dynamics: Dashed and dotted lines in Fig. 6 show corresponding thermalization dynamics at T = 30 K and T = 77 K, respectively, exhibiting population decay. The temperature dependence clearly shows that correlation lengths within the full system-reservoir dynamics are of importance, and the observed effect allows to probe these otherwise inaccessible microscopic environmental properties. This formation of self-stabilizing dissipative structures is closely related to a localized phase stabilization in the coherent driving case Ω 0 = 0 at Ohmic spectral densities and without photons [31,43]. There, above a critical coupling strength the system transitions into a localized phase with nonzero steady state population rather than decaying to zero. In our case, a similiar phenomenon is established with incoherent feedback instead of coherent external driving, addressing the localized phase stabilization process from a dissipative non-Markovian side. The complexity of this phenomenon is illustrated in the inset in Fig. 6, showing the dimension of the link index connecting the current open system bin to the discrete memory MPS over time (see Fig. 4). The red line shows the case including the continuous reservoir. After slowly increasing during times t < τ , it exhibits a vast increase once feedback sets in and quickly reaches a maximum due to finite memory. The high grade of entanglement between the two reservoirs even at arguably low finite memory sizes n d = 4, n c = 4 underlines the crucial role of the interplay between the two non-Markovian processes for the observed protection of coherence. Switching off the phonon coupling, g q = 0, results in a much lower maximum link dimension (blue line), as no entanglement between the reservoirs arises. For phenomenological dephasing (orange line), the link connecting system and memory bins has an even lower dimension due to the highly decreased complexity of the then 1D network.

V. CONCLUSIONS
We have presented an MPS algorithm for the description of a time-discrete quantum memory in Liouville space. By combining this technique with a path integral tensor network implementation for continuous non-interacting harmonic reservoirs, we have established a quasi-2D tensor network, allowing for simulations of quantum systems subject to two non-Markovian environments while maintaining crucial entanglement information in the coupled system with both diagonal and off-diagonal system-reservoir interactions. Due to arising reservoir-reservoir correlations, system correlations scale intensively with respect to the twofold single reservoir case. In consequence, the achievable memory depth is limited to n c + n d < 20 in our study. However, appropriate tuning of the relevant system and reservoir time scales via the employed parameters still opens up a wide array of accessible systems and scenarios where numerical convergence can be achieved. The next step will be to trace out the time-discrete feedback bins as well after their interaction with the systems' degrees of freedom to further improve numerical efficiency. This will allow for longer delay times and therefore increased time discretizations without changing the qualitative results. Hence, the presented quasi-2D tensor architecture is a first step towards  Fig. 6 with respect to (a) the continuous reservoir memory depth nc and (b) the Schmidt value cutoff precision dcut applied during the singular value decomposition. In (a), the left inset shows a zoom-in on the long term dynamics, with the right inset depicting the resulting dimension of the link index connecting the two tensor networks. The inset in (b) shows a zoom-in on the long term dynamics.
unraveling the mostly unexplored field of multiple interacting non-Markovian reservoirs in a numerically exact fashion.
As an example application, we have demonstrated that the interplay of a structured phonon reservoir and photon feedback can dynamically protect the system from destructive interference by time-delayed backflow of coherence, resulting in dynamical population trapping. Tuning the non-Markovian interactions with respect to each other allows for the formation of inter-reservoir correlations, dynamically preservering feedback-induced coherence in the system. These findings have implications for the fields of quantum thermodynamics and nonequilibrium physics, as well as dynamical quantum phase transitions with ergodic, entropic or negentropic information exchange, where taking account of such dissipative structures may unravel new phenomena. Future works will aim to advance our architecture to a full 2D representation via projected entangled pair states with combined memory bins, potentially allowing for simulations of multi-level systems at improved time resolutions.  Fig. 2) and (c) the combined scenario (see Fig. 6). (d) shows a convergence analysis with respect to the order of the series expansion performed in the Liouvillian L(t), with the inset showing a zoomin on the long term dynamics.

Appendix A: Convergence analysis
Here we provide a detailed analysis on the numerical convergence of the presented results. In light of the memory limitations imposed by an intensive scaling of entanglement when describing two non-Markovian reservoirs simultaneously, we first confirm convergence with respect to the memory depth of the continuous harmonic reservoir n c , i.e., the validity of the finite memory approximation for the considered scenario. Fig. 7(a) shows the population dynamics corresponding to Fig. 6 at ϕ = 1.17, T = 4 K and increasing memory depths n c of the contin-uous reservoir. The left inset shows a zoom-in on the long term dynamics, while the right inset depicts the resulting dimension of the link index connecting the two tensor networks, demonstrating the here occurring exponential growth of inter-reservoir entanglement with increasing memory depths. While the results are clearly not convergent at n c = 2 (dark blue line), the cases n c = 3 and n c = 4 (light blue and green lines) already are in very good agreement. Between the cases n c = 4 and n c = 5 (dashed orange line), no difference can be seen even at close range (see left inset), hence justifying the choice of n c = 4 employed in our calculations. Fig. 7(b) shows the same dynamics calculated at increasing Schmidt value cutoff precisions d cut applied during the singular value decomposition [46]. Here, no differences between the cutoffs d cut = 10 −8 (dark blue line) and d cut = 10 −14 (dashed orange line) can be seen even at close range (see inset), underlining the convergence of our results with respect to the employed cutoff precision d cut = 10 −12 (green line). Hence, no crucial entanglement information has been truncated during the time evolution.
In addition, we investigate the numerical convergence of our results with respect to the time evolution step size ∆t to ensure that no errors occur during the Trotter decomposition. We first calculate the dynamics of both reservoirs independently for decreasing step sizes. Fig. 8(a) shows the time evolution dynamics resulting from the time-discrete reservoir at ϕ = 1.17 for decreasing step sizes ∆t (see Fig. 5). Aside from minor differences during the initial time steps, the long term dynamics show very good agreement for all employed ∆t, underlining the convergence of the feedback algorithm at the chosen step size ∆t = 0.3 ps (light blue line). Fig. 8(b) shows the polarization dynamics imposed by the continuous reservoir [see Fig. 2(b)], calculated at T = 4 K and decreasing time evolution step sizes ∆t. Again, the resulting long term dynamics at ∆t = {0.3 ps, 0.15 ps, 0.075 ps} are in good agreement. Since both processes converge individually with respect to the time discretization, the combined setup can be expected to converge as well, since the involved time scales remain the same. Fig. 8(c) shows the dynamics of the combined system coupled to both structured reservoirs for time discretizations ∆t = {0.3 ps, 0.24 ps}, corresponding to n d = {4, 5} feedback bins. Due to limited computational resources imposing restrictions on the memory depths of both reservoirs and thus the maximum achievable time discretization at a given feedback time τ , the presented calculations are performed until t = 10 ps, where they exhibit good agreement with each other. Lastly, in Fig. 8(d) we verify the convergence of the numerical implementation with respect to the order of the series expansion performed for the Liouvillian L(t) [see Eq. (13)]. While the resulting dynamics show minor variations between the eighth and ninth order series expansions on a close scale (blue and green lines, see inset), no difference can be observed when comparing the dynamics resulting from ninth and tenth order expansions (dashed orange line). In conclusion, the employed tenth order series expansion of L(t) ensures convergent results as well.