Information Loss Pathways in a Numerically Exact Simulation of a non-Markovian Open Quantum System

In the non-Markovian regime, the bath has the memory about the past behavior of the open quantum system. This memory has slowly-decaying power-law tails. Such a long-range character of the memory complicates the description of the resulting real-time dynamics on large time scales. In a numerical simulation, this problem manifests itself in a `revival', a spurious reflected signal which appears after a finite time thus invalidating the simulation. In the present work, we approach this problem and develop a numerical discretization of the bath without revivals. We find that a crucial role is played by the singularities of the bath spectral density (e.g. edges of bands): the memory about the (spectral) behaviour of the open system in the remote past is completely averaged out (forgotten), except an increasingly small vicinity of these singular frequencies. Therefore, we introduce the concept of memory channel, to denote such an irreversible information loss process around a particular singular frequency. On a technical side, we begin the treatment by noting that with respect to the memory loss, the quantum field of the bath should be divided into the following two different parts: the observable and the virtual quanta. Information about the former is never lost: they contribute to the trace over the bath degrees of freedom. The only way to avoid the revival from the observable quanta is to calculate their dynamics exactly. We do this by employing a stochastic sampling of the Husimi function of the bath. The other part, the virtual quanta, are always annihilated after a certain delay time. We construct a dedicated quantum representation for the virtual states by assigning amplitudes to the delay times before annihilation. It is in this representation we rigorously define the concept of memory channels.

In the non-Markovian regime, the bath has the memory about the past behaviour of the open quantum system. This memory has slowly-decaying power-law tails. Such a long-range character of the memory complicates the description of the resulting real-time dynamics on large time scales. In a numerical simulation, this problem manifests itself in a "revival", a spurious reflected signal which appears after a finite time thus invalidating the simulation. In the present work, we approach this problem and develop a numerical discretization of the bath without revivals. We find that a crucial role is played by the singularities of the bath spectral density (e.g. edges of bands): the memory about the (spectral) behaviour of the open system in the remote past is completely averaged out (forgotten), except an increasingly small vicinity of these singular frequencies. Therefore, we introduce the concept of memory channel, to denote such an irreversible information loss proccess around a particular singular frequency. On a technical side, we begin the treatment by noting that with respect to the memory loss, the quantum field of the bath should be divided into the following two different parts: the observable and the virtual quanta. Information about the former is never lost: they contribute to the trace over the bath degrees of freedom. The only way to avoid the revival from the observable quanta is to calculate their dynamics exactly. We do this by employing a stochastic sampling of the Husimi function of the bath. The other part, the virtual quanta, are always annihilated after a certain delay time. We construct a dedicated quantum representation for the virtual states by assinging amplitudes to the delay times before annihilation. It is in this representation we rigorously define the concept of memory channels. The description of a large number of physical phenomena is based on the picture of a bounded quantum subsystem (the open quantum system, OQS) which is coupled to a continuum (the bath). These include the divesre physical scenarios like the atoms or quantum dots coupled to a continuum of modes of phononic, electronic, and photonic reservoirs [1]; resonance phenomena in the electron-molecule collisions [2][3][4][5], where the colliding electron and the molecule form a joint interacting intermediary state (the bounded subsystem), and the continuum is formed by outgoing (scattered) states. In the reactive molecular collisions [6], the intermediary state is coupled to the two continua: the incoming states of reactants, and the outgoing states of products. In all the cases the total problem has infinite dimensions, and due to the interactions present in the subsystem, it has no trivial solution. The traditional approach is to place a sharp boundary between the subsystem and the continuum: the subsystem is interacting and nontrivial, but finite; the continuum is infinite, but simple. Then, by (approximately) eliminating the continual degress of freedom, one obtains a reduced finite description, which is amenable to numerical solution. This is the essence of the Feshbach projections method [3,[7][8][9].
There is no problems with the Feshbach projections in the Markovian mode. This regime is characterized by the two conditions: 1) the coupling between the subsystem and the continuum is weak, and 2) the continuum is much faster then the subsystem. The condition 2) is also may be stated in a different way: that the background (continuum) spectral density is flat on the scales of spectral features of the subsystem [1]. In operational terms, this often means that we need to find such a boundary (between the subsystem and the continuum) which fullfils these conditions [2]. In the case of a success, the boundary between the continuum and the subsystem is crossed instantly and only once by the incoming/outgoing particles (bath excitations). This results in a simple reduced description. For the open quantum systems, like atoms and quantum dots, the reduced description amounts to the addition of a Lindbladian dissipator to the von Neumann equation of the reduced density matrix [1]. In the case of resonant and reactive collisions, a local optical complex potential [2,6] is introduced into the Hamiltonian. Thus, nowadays we have a fairly complete and consistent picture of the Markovian mode.
However, outside the validity of this mode, in the non-Markovian regime, the concept of sharp boundary becomes problematic. The continuum acquires a memory about the past subsystem behaviour. In other words, there is non-zero amplitude that the particles (bath excitations) will return back from the contuum to the subsystem: the emitted virtual photons can be reabsorbed by the atom [1]; the reaction fragments can temporarily return to the interacting region during a reative collision [10][11][12][13][14]. Moreover, the memory of the continuum is long-range: no matter how far the particles fly into the continuum, the return amplitude decays only as inverse power law. This means that physically there is no sharp boundary at all. And when we artificially place a sharp boundary, the low energy spectrum of the problem will be distorted [2][3][4]. In the time-resolved simulations, this corresponds to the revivals after a finite time [15][16][17], so that large-time asymptotical behaviour of the observables is corrupted. When doing numerical computations, one may attempt to alleviate this problem by shifting the boundary inside the continuum e. g. either by a direct inclusion of the continuum modes or sites, or by shifting the memory tail cutoff in the time-nonlocal approaches (quasi adiabatic path integrals and related [16][17][18]). But then, the dimension of the reduced description becomes combinatorially large, and we cannot advance much this way.
In summary, the long-range memory tails are the major obstacle for the efficient description of low-energy and large-time physics of open quantum systems. In the present paper we propose one possible solution to this problem.
Our treatment is based on our preceding paper [19] where we have presented an alternative approach to the description of open quantuam systems in a non-Markovian bath: we do not place a sharp boundary between the retained and the eliminated degrees of freedom in a reduced description. Instead, the quantum field of bath excitations is divided into the following two components: 1) the excitations which are irreversibly (only once) emitted by OQS, and which are observable by a bath measurement, and 2) the virtual excitations which are always reabsorbed by OQS before any measurement of the bath [20,21]. The former component, the observable quantum field, has a classical stochastic structure, since its time evolution can be described by a master equation for the probability distribution of the measurement outcomes of the bath. Therefore, we can numerically simulate it by Monte Carlo methods. The latter component, consisting of virtual excitations, is a genuinely quantum object. Its time evolution should be computed by a deterministic scheme. In other words, we obtain a reduced description in which the virtual part of the bath state is retained, whereas the observable part of the bath state is eliminated by addition of a classical stochastic noise. One can understand this as a "soft boundary" between the OQS and the bath.
The two components of the bath quantum field differ from each other in how they carry the memory about the past OQS behaviour. The observable quanta directly contribute to the measurement, thus to the reduced density matrix of OQS. Therefore, the information in the observable quanta is not lost. In order to avoid the revivals from the observable component, we propose to calculate its evolution exactly, by a Monte Carlo simulation. At the same time, in the second component, the virtual quanta, only the information in the return amplitudes is retained. Therefore, there is a mechanism of the gradual loss of memory about remote past. We identify this mechanism: it is based on a monotonically increasing flattering of the memory tails [16]. This allows us to devise a "soft" discretization of the bath, so that all the effects of the long-range tails are taken (numerically) into account.
Our approach is supported by the current developments in the literature on simulation methods for the timeresolved simulations. In particular, there are numerous combined stochastic-deterministic schemes [22,23]. In these algorithms, the interaction of the bath with OQS is splitted into the two parts: one is represented by a stochastic coloured noise (whose correlator reproduces the bath memory function); the other part is solved deterministically. These algorithms have demonstrated promising convergence properties when simulating non-Markovian OQS, at short to moderate time scales [22,23].
In section II A we introduce the model of open quantum system. We define the problem of long-range memory tails in section II B. Then, in section II C we briefly recapitulate our approach of observable and virtual quanta. We demonstrate that the observable quanta can be stochastically simulated without the revivals in section II E. Then, in section III we develop a special representation for the virtual states which helps us to identify the memory loss mechanism in section III C. In particular, we find that the memory loss is concentrated around singularities of the bath spectral densities, which we call "memory channels". This allows us to implement a "soft" discretization of the bath, whose convergence properties are illustrated on a calculation in section III E We conclude in section IV. In appendix A we derive the equation of motions for OQS in the dedicated representation for virtual states. In appendix B we provide details of how the channelwise memory functions are constructed. Finally, in appendix C we provide some detials about our implementation of the numerical methods

II. MEMORY OF THE QUANTUM ENVIRONMENT
We begin the exposition by introducing the model of open quantum system, formulating the problem of long-range memory tails, and recapitulate the approach of dressed quantum trajectories to it.

We consider a model of open system which is bilinearly coupled to the bath of harmonic oscillators. The Hamiltonian is
where H s is the OQS, H b is the bath The coupling is through the operator s which is in the system's Hilbert space, and through the operator b which is in the Hilbert space of bath,

B. Long-range memory tails
The memory function of the bath is defined in terms of the commutator of the coupling operator (in the interaction picture) For physical bath couplings (with no negative frequencies) the memory function always has inverse power-law tails. This can be illustrated by considering a typical form of the spectral coupling coefficient where α > 0 is the coupling strength; for s < 1 we have the so-called subohmic bath; s = 1 corresponds to ohmic bath; finally, the case s > 1 is called the superohmic bath. The corresponding memory function is C. Two types of bath quanta Our goal is to devise a simulation scheme which takes into account the full memory tails in a numerically exact way. Such a scheme would enable us to perform long-time simulations without the revivals. Starting from an initial factorized state the open system begins to interact with the bath. This interaction manifests itself in an emisstion and absorbtion of quanta (bath excitations). Then, after a time t, we perform a non-selective measurement of the bath state, which yields the reduced density matrix for the open system. Here the partial trace operation Tr b is done with respect to the bath degrees of freedom. From the point of view of non-selective measurement (the Tr b ), all the bath quanta emitted by OQS are divided into the following two types: those which will survive up to the measurement time t (the observable quanta), and those which will be absorbed by OQS (the virtual quanta).

D. Exact simulation of observable quanta
The memory about the observable quanta is not lost, since these quanta contribute to the Tr b and to the reduced density matrix. There is no obvious way to take into account this "observable" memory on long times except a certain numerically exact simulation. Any approximate simulation (e.g. finite discretization of the bath) will yield strong revival signals. Due to infinite dimension of the bath Hilbert space, a Monte Carlo simulation is the most appropriate numerically exact approach. As was discussed in the preceding paper [19], the natural way to construct such a stochastic simulation is to probabilistically sample the outcomes of a bath measurement. By considering the measurement with respect to the coherent states of the bath, in [19] the following Monte Carlo simulation algorithm was presented. At a time moment t, the reduced density matrix ρ s (t) of OQS is represented as a statistical average of pure states, where there averaging is over the complex Gaussian spectral noise z (ω), with the the statistics For each realization of the noize z (ω), the pure state |Ψ dress (z, t) is a solution of the Schrodinger equation where the noise-dependent Hamiltonian is Here, the time-dependent classical complex signal ξ (t) depends on the noise,  Figure 1. In the previous paper [19], a simulation of the driven two-level system in a waveguide was performed. n = 0 curve corresponds to the case of zero virtual quanta: no revivals at extremely large times. n = 1 curve corresponds to the case of inclusion of one virtual quantum: a spurious revival signal is present at times after t = 400.
Another complex time-dependent classical field self-consistently depends on the average of the coupling operator s at previous times.

E. Problem of memory tails for the virtual quanta
Such an approach -the exact Monte-Carlo simulation of the observable quantum fields -solves the revival problem for these observable fields. In order to demonstrate this, in our previous paper [19] we performed the numerical simulation of the driven two-level system in a waveguide.
The Hamiltonian for the driven OQS is with the driving field f (t) = 0.1 cos t. The spin is coupled to the waveguide through the spin lowering operator s = σ − . The degrees of freedom of the waveguide are the spectral modes labeled by wavevector k ∈ [0, π], with frequencies ω (k) = ε 0 + 2h cos k and the coupling constant c (k) = h 2 π sin k. In Fig. 1 we present the simulation results for h = 0.05, ε = 1 on large time scale. In this figure, n = 0 curve corresponds to the approximation when all the virtual quanta are neglected (discarded all the operators b (t), b † (t)) in Eq. (12). We see that in this case, when there is only the observable quantum field, the system reaches the steady state, and no revivals are present.
At the same time, if we impove this approximation and include one virtual quantum in the simulation, which corresponds to n = 1 curve in Fig. 1, and simulate the one-quantum part of the dressed wavefuction Ψ dress (z, t) on a discretized grid of N = 20 modes, we observe a clear revival signal. See [19] for details.
We conclude that the stochastic simulation of observable field solves only one half of the memory tail problem. The other half, for virtual quanta, remains. In the next section we propose a solution to this problem.

III. AMPLITUDES FOR VIRTUAL QUANTA
In this section we deal with the question of how to efficiently compute the projection at long times. Usually, the state |Ψ dress (z, t) is represented as a vector of probability amplitudes over a certain bath basis. The conventional choices for the basis are either the frequency modes, or the sites of the equivalent semiinfinite chain [24][25][26]. However, the disadvantage of this approach is that whenever we truncate the basis, a spurious reflected signal (the revival) will come, thus invalidating the simulation after a finite time. In order to delay this revival (in other words, to increase the simulation time), the number of frequency modes (or sites) should be increased linearly with time. This makes the long-time simulation prohibitive, especially if a rather large number of virtual quanta is retained. The solution to this problem begins with the following observation. It is inefficient to perform the computation exactly according to (17). Because this way we first solve for the full state |Ψ dress (z, t) , which contains both the observable and the virtual quanta. And only when the observables are computed, we project it to 0 b |, thus discarding the observable quanta. Therefore, instead of |Ψ dress (z, t) , it is reasonable to find such a representation which contains only virtual quanta from the outset, and derive the equations of motion directly in this representation.

A. The delay-time amplitudes
The characteristic property of the virtual quanta is that each of them will be absorbed after a certian delay time. Then, assume we have the bath in a superposition of number states For each number of quanta N , we introduce a set of delay times τ 1 , . . . τ N , and define the delay-time amplitudes as . . .
Here in ϕ N (τ 1 , . . . , τ N ; t) we use the semicolon in order to distinguish the external (laboratory) time t from the dynamical variables -the delay times τ k . The amplitude ϕ N (τ 1 , . . . , τ N ; t) is a symmetric function of its arguments.

B. Equations of motion in the delay-time picture
Now let us write the equation for the dressed quantum trajectory (11)- (15) in terms of the delay-time amplitudes. In this case, besides the bath degrees of freedom (now represented by delay times), we also have the OQS quantum numbers. Denoting these additional quantum numbers as s, we have the joint amplitude ϕ N (s, τ 1 , . . . τ N ; t) that all the virtual quanta will be absorbed after the delay times τ 1 , . . . , τ N , and OQS will be in the state |s : We can also look on this from the other side, that we have a ket vector |ϕ N (τ 1 , . . . τ N ; t) s in the OQS Hilbert space, which depends on delay times, and which is related to the joint amplitude as In order to derive the equations of motion for these delay-time ket vectors, we differentiate with respect to time their definition (23): Evaluating different terms in H dress , we obtain a hierarchy of equations, which we present here up to the second order where we have grouped into H stoch all the terms which belong to the Hilbert space of OQS: The average values of s are now computed as We refer the interested reader to the appendix A for the details of the derivation.

C. Irreversibile loss of memory at large delay times
Now let us analyze the equations (26)-(28). All the information which is not lost but which is accessible to the outside world, is contained in the shifted classical noise ξ (t) + φ * (t). At the same time, the virtual quanta is a mechanism of the gradual loss of information. In order to see how is it working, let us have a look at a typical memory function, Eqs. (5)-(6). In Fig. 2 we plot real part of M (τ ), in logarithmic time scale. Imaginary part has qualitativly similar behaviour. What we notice is that at large delay times the tails of the memory function become increasingly flat. More formally, we can say that with the increasing delay τ , there is an increasingly large time scale ∆τ such that the memory function can be considered effectively constant on the intervals [τ, τ + ∆τ ]. However, since in the equations (26)-(28) the system states s |φ 0 (t) and s |φ 1 (τ i ; t) are continuously superimposed in the bath delay variables, with the non-local weight M (τ ), this means that at large delays the information about OQS trajectory on the time scales smaller then ∆τ , is averaged out (irreversibly lost). And asympotically, at infinitely large delays, only the information about zero-frequency behaviour of the system is retained. Therefore, the long-range irreversible memory-loss mechanism is based on the monotonically increasing flattering of the memory tails.

D. Case of a finite bandwidth
Now let us return to the system treated in the section II E. The semiinfinite bath for that model has the memory function The behaviour of the tails in this case is a little bit harder. From Fig. 3 it is seen that while descreasing, the memory tails always oscillate, and we cannot apply here directly the reasoning from the previous section: there is no flattering of the tails. However, if we review other numerical methods which employ the ideas of coarse-graining (e.g. the numerical renormalization group method of Bulla et al. [27]), we notice that in these methods all the memory  (long-time behaviour of the bath) is concentrated around the edges of the band, which is reflected by the fact that we need exponentially-fine descrerization of the bath spectrum at these points. And the oscillations in Fig. 3 are just the manifestation of the interference between these two memory channels. In other words, we expect that if we write the memory function separately for each channel, then it will again have simple, increasingly-flattering tails. Formally, we can define the notion of memory channels through the large-delay asympotical behaviour of the memory function. In the case of Eq. (31) we have Therefore, indeed, we have the two memory channels, associated to the two momotonically-flattering power-law tails. Each channel is attached to the frequency of the corresponding band edge, ε 1 = ε − 2h, and ε 2 = ε + 2h.
We expect that this feature is general: each edge of the band, and each singularity in the spectral density, will lead to the existence of its own channel. In order to make explicit this memory-channel structure, we decompose the memory function: where and R (z) is a certain fast-decaying real-valued function, which fixes τ = 0 singularity of Hankel functions. See Appendix B for the details about R (z). Let us rewrite the equations (26)-(28) in terms of the channel memory functions. Since we now have the two memory functions, the delay-time amplitudes get additional indices, indicicating to which channel the delay argument referes to. It can be interpreted as an additional "quantum number" of the virtual quantum. So we obtain the set of wavefunctions: ϕ 0 (t); φ For one virtual quantum in the kth channel we get Finally, the equations for the two quanta in the same channel k, and for the two quanta in different channels k = l The numerical coarsegraining of delay-time amplitudes Now its time to implement our understanding of the irreversible memory loss mechanism. Since the time-scales of the memory functions increase at large delay-times, we "compress" them by making a substitution where in our numerical calculations we employ a = 0.5, b = 0.1, c = 0.9. Actually, we didn't try to find the most optimal substitution. The form (41) was guessed by eye: so that the memory function is not too compressed near τ = 0, and at the same time, that the power-law tails are sufficiently pursed. Then, the memory function is considered as a function of the independent argument ς: To express the equations (37)-(40) in terms of the ς-variables, we make the following substitutions: and A grid of ς values was employed: where and m = 7.
We have employed the cutoff value All the amplitudes depending on delay times were approximated as piecewise polynomials on the intervals [ς i , ς i+1 ].
The order and the coefficients of the polynomials where determined from the requirement to reproduce the values and first three derivatives of the amplitudes on the grid ς i . Yet, we don not argue that this is the most efficient coarse-graining and interpolation scheme. Rather, in this work we provide it as a proof-of-principle.
The resulting system of equations was solved by a second-order symmetric split-operator approach, where we considered separately the propagations generated by free shift in delay time, by hopping with the memory function M k (τ ), by absorbtion of virtual quantum, and by the Schrodinger evolution. See the Appendix C for the details about the slip-operator method. The results are presented in Fig. 4 From this graph we see that the system successfully reaches the stationary non-equilibrium state. There is no revivals or other deterioration of the accuracy of observables.

IV. CONCULSION
In this work we present a solution for the problem of how to discretize the bath so that there is no revivals at large times. We have found a soft coarse graining procedure for the memory function so that the whole effect of its long-range tails is taken into account, in a numerically exact way. We illustrate this approach by providing the results of test calculations for the driven spin-boson model: only two virtual excitations are enough to achieve the uniform physical convergence on a large time interval.
Another interesting result is that each singularity in the spectral density of the environment (e.g. band boundary) leads to the existence of its own memory channel, which "operates" on the frequency of the singularity. The physical content of this result is that at large delay times, the memory about the past behaviour of OQS is retained only in the vanishingly small vicinity of these singular frequencies, and is completely averaged out (forgotten) away from this discrete set of spectral lines. Now the actual task is to extend this result to a strong coupling regime, where a large number of virtual quanta may be present. One of the prossible approaches is to implement a singular value decomposition of the coarse-grained delay-time wavefunction of virtual quanta, in order to further "compress" the virtual quanta wavefunction. This is a matter of future investigation.

ACKNOWLEDGMENTS
The study was founded by the RSF, grant 16-42-01057. . Simulation of the driven two-level system in a wave-guide, described in section II E. Two virtual quanta are taken into account. When a soft coarse-graining of the memory function on large delay times is performed, the revivals dissappear.
There is no deterioration of accuracy on large time scales.

Appendix A: DERIVATION OF EQUATIONS IN THE DELAY-TIME PICTURE
In order to derive the equations of motion for the delay time amplitudes, we need to consider the following three ingredients of H dress (z, t): the free bath motion, the action of the annihilation operator, and the action of the creation operator.
Let us derive the equation of motion for the delay-time amplitudes under a free evolution of the bath: Here, in the first line we have employed the Schrodinger equation for the free bath state |ϕ N (t) . Then, we moved H b to the left by commuting it through all b (τ p ). Finally, in order to arrive at the last line, we use the properties To describe the action of the annihilation operator b, we introduce a "quantum-subtracted" bath state In order to derive the representation of the last element, the action of the creation operator, we introduce a "quantumadded" bath state and employ the two-time canonical commutation relation (4): . . . Now, expanding the definition of H dress (z, t) in the equation (25), and substituing the derived here relations, one obtains the result (26)-(28).