A time-dependent regularization of the Redfield equation

We introduce a new regularization of the Redfield equation based on a replacement of the Kossakowski matrix with its closest positive semidefinite neighbor. Unlike most of the existing approaches, this procedure is capable of retaining the time dependence of the Kossakowski matrix, leading to a completely positive divisible quantum process. Using the dynamics of an exactly-solvable three-level open system as a reference, we show that our approach performs better during the transient evolution, if compared to other approaches like the partial secular master equation or the universal Lindblad equation. To make the comparison between different regularization schemes independent from the initial state


Introduction
Describing the time evolution of a quantum system interacting with an external environment is of paramount relevance in contemporary physics, with applications in a wide variety of fields.However, portraying the exact dynamics of an open system is often challenging, due to the intrinsic complexity of dealing with the large number of degrees of freedom of the environment.This problem can be tackled by introducing effective descriptions in which only a small, but essential, amount of bath properties is taken into account to derive a good approximated picture of the system evolution.One of such examples are Markovian Quantum Master Equations (QMEs), which express the time derivative of the density matrix in terms of a superoperator satisfying the prescriptions of the Lindblad-Gorini-Kossakowski-Sudarshan (LGKS) theorem [1,2].
QMEs are one of the most widely used models in open quantum systems and have been applied to problems of quantum transport, computation, chemical modeling and quantum thermodynamics [3,4].Despite that, the range of validity of QMEs and their reliability in the description of coherent effects, are still debated today [5].The standard Born-Markov approximation [6,7] leads to the Redfield equation [8,9], which is known to violate the positivity of the density operator [10] and other desirable properties of the quantum evolution [11,12], hence providing nonphysical predictions.Historically, the main route for curing these issues required the energy levels of the system to be well separated: this is the secular approximation, which provided remarkable results in the context of quantum optics and quantum chemistry [13,14].This approach is however not suited for the study of generic many-body systems, where the spacing between energy levels typically decreases exponentially when increasing the size.For this reason, in recent years a number of works appeared in the literature proposing ways to obtain LGKS equations that are free of the restrictions imposed by the secular approximation [15][16][17][18][19][20][21][22][23][24][25][26].
In this work we argue that these "regularization" techniques amount to a substitution of a certain matrix that describes the dissipative dynamics, known as the Kossakowski matrix, with a positive semidefinite one, thus leading to a LGKS equation.With this observation at hand we propose a natural regularization scheme for the Redfield equation, which consists in replacing the Kossakowski matrix with its closest positive semidefinite matrix of the same dimension.The result is compared with some of the existing schemes, by examining to what extent they reproduce the true dynamics of a simple open system that can be solved exactly.To do this, we employ a novel technique that makes use of the Choi-Jamiołkowski isomorphism [27,28] to envision a numerical comparison that is independent of the choice of the initial state of the evolution.We emphasize that our procedure can be applied not only to the standard Redfield equation but also to its version with time-dependent coefficients (that can result, for instance, from avoiding the so-called "second Markov" approximation).In this case, our regularization preserves the time dependence of the coefficients but makes the dynamics completely positive divisible [29][30][31].The residual time dependence is an indicator that our approach could perform better at short times with respect to existing schemes.
The paper is structured as follows.In Sec. 2 we write the Redfield equation and the associated Kossakowski matrix.In Sec. 3 we discuss the regularization procedure, first by analyzing common existing schemes in Sec.3.1 and then by presenting our natural proposal in Sec.3.2.In Sec. 4 the reader can find an example of application to an open three-level system: in Sec.4.1 we present a direct numerical calculation at the level of the density matrix, while in Sec.4.2 we discuss the novel Choi operator technique to compare predictions of different master equations.Finally, in Sec. 5 we draw our conclusions.

The LGKS theorem and the Redfield equation
Let us consider a quantum system S described by a Hilbert space H S of finite dimension N .Moreover, let L(H S ) be the vector space of linear operators on H S equipped with the Hilbert-Schmidt inner product 〈X , Y 〉 := Tr X † Y .A linear map Φ : L(H S ) → L(H S ) is expected to describe a physical transformation when it is trace-preserving and completely positive (CPT) [6].
A particularly important case is when the process is described by a quantum dynamical semigroup [32], i.e., a one-parameter family {Φ t } t≥0 of linear CPT maps on L(H S ) such that t → Φ t is continuous, Φ 0 = 1, and Φ t+s = Φ t • Φ s .In this case, if ρ(0) is the initial state of S then the state at time t is given by ρ(t) = Φ t (ρ(0)).Quantum dynamical semigroups are important because given {Φ t } we can find a linear operator L, called the generator of the semigroup, such that dρ(t which is in the form of a Markovian master equation [6].The well-known theorem of Lindblad [1], Gorini, Kossakowski, and Sudarshan [2] characterizes the shape of such a generator.Here we will use the formulation provided by [33], which is best suited for our discussion. Theorem.A linear operator L : L(H S ) → L(H S ) is the generator of a quantum dynamical semigroup if it can be written in the form where H = H † , {F i } i=1,...,N 2 is an orthonormal basis of L(H S ) and χ is a positive semidefinite complex matrix which is uniquely determined by the choice of {F i }, called Kossakowski matrix.
One could be also interested in generalizations of Eq. ( 1) in which the generator becomes time dependent L t .In this kind of scenario we must deal with a two-parameter semigroup Φ t,s = T exp t s dτ L τ , where T is the time ordering.Here divisibility is guaranteed, in the sense that for any t ≥ s ≥ 0 there exists a linear map Λ t,s (intertwining map) such that Φ t,0 = Λ t,s • Φ s,0 .The case of completely positive intertwining maps (CP divisibility) is particularly important, since it is characterized by a LGKS-like generator (2) where χ and H become time-dependent quantities [29][30][31].At this level, notice that χ(t) ≥ 0 for all t ≥ 0 is a sufficient condition for CPT dynamics but it is by no means necessary [15].
In a general setting where S is allowed to interact with an environment E, the non-Hamiltonian terms in Eq. ( 2) play a crucial role.This can be seen with a microscopic derivation of (2), in which we start from a unitary description of the universe U = S ∪ E, trace away the environment and obtain, under suitable assumptions, a master equation for S which is in LGKS form.Since the universe is closed by definition, it is described by a Hamiltonian, which is commonly written as The interaction term H I is of the form where A α acts on the system and B α acts on the environment.It is also common to assume these coupling operators to be Hermitian, but here we will not make this assumption.
Under the Born-Markov approximation (and other standard assumptions) it is known that the reduced dynamics of S in the interaction picture ρ(t) := e iH S t ρ(t)e −iH S t is [6] where c αβ (τ) := 〈 B † α (τ)B β 〉 is the environment correlation function (the average 〈•〉 is calculated on the stationary state of the environment).This is one of the many forms of the so-called Redfield equation [8,9], and it is obtained under the assumption that the typical evolution time τ S of ρ(t) is the longest timescale of the problem (see [20] for further discussions).
Our first task is to write the Redfield equation in the form (2). Let us consider the basis {|k〉} of normalized eigenvectors of the free system Hamiltonian H S , so that we can write its spectral decomposition as H S = k ω k E kk , where E kq := |k〉〈q|.Since {E kq } k,q=1,...,N is an orthonormal basis of L(H S ) we can expand, for example, where A β,kq = 〈k|A β |q〉.We have now to replace the decomposition (6) inside Eq. ( 5): the details are reported in App. A. Going back to the Schrödinger picture, one finds where and We also introduced the Bohr frequencies ω kq := ω q − ω k .Since χ and H LS generally depend on time, we do not have an equation of the form (2) yet.A common way to strictly obtain an equation in the form (2) is to perform the second Markov approximation [6], which consists in replacing in (10).In this scenario, we will write Γ αβ (ω) = lim t→∞ Γ αβ (ω, t).For future reference, notice that if we split Γ αβ = J αβ + iS αβ in its real and imaginary part, we have where ĉαβ is the Fourier transform of c αβ .Therefore one can invoke Bochner's theorem to infer that J is a positive semidefinite matrix [6].
In the general time-dependent case it is difficult to say when the Redfield equation leads to CPT dynamics, since (to our knowledge) there have not yet been found necessary conditions for complete positivity of a time-dependent generator.However, we do know a sufficient condition, namely χ(t) ≥ 0 for all t ≥ 0. Unfortunately (as we shall also see in Sec.3.2), χ(t) is not positive semidefinite in general for the Redfield equation.In the time-independent case, this is a long-standing well-known problem [10].

Common existing regularizations
We argue that most of the procedures to recover positivity from the Redfield equation are in fact ways to transform the matrix χ in Eq. ( 8) into a positive semidefinite one, and the vast majority of them only deals with the time-independent case.For example, a popular approach consists in performing a coarse-graining transformation in the equation for ρ(t) [18].In our notations, this means that we apply the operation to both sides of Eq. ( 45).The advantage lies in the fact that if we choose ∆t ≪ τ S , where τ S is the timescale of variation of ρ(t), we can ignore the action of C ∆t on ρ(t) and pull the latter out of the integral.Using the fact that where sinc(x) = sin(x)/x is the cardinal sinus, we see that the effect of the coarse graining is given by the substitution and similarly for the Lamb shift coefficient η kq,nm → η kq,nm .The interesting fact about this expression is that one can prove that if ∆t is sufficiently high the matrix χ (∆t) will be positive semidefinite [18].In the extreme situation ∆t → ∞ one has which is the Kossakowski matrix obtained with a secular approximation [6].For this reason one also says that χ (∆t) with general (but appropriate) ∆t is the Kossakowski matrix in partial secular approximation.
A more recent and permissive construction is the one provided by Nathan and Rudner in Ref. [22] and by Davidović in Ref. [5].The idea is to replace the arithmetic mean that appears in ( 8) with a geometric one: where it is intended matrix multiplication of matrix square roots (recall that J ≥ 0).It can be shown that this approach is justified whenever τ S ≫ 1/ω R , where ω R is representative of the frequency range of the system [34].

New regularization
Since we have to regularize the Kossakowski matrix, the following proposal seems very natural: for every t ≥ 0, replace χ(t) with its closest positive semidefinite matrix of the same dimension.More precisely, given a norm ∥•∥ on the space of N × N complex matrices we define χ + (t) := arg min and use χ + (t) instead of χ(t) in the Redfield equation, thus obtaining a LGKS-like equation.Unlike standard approaches, notice that here we are retaining the time dependence of χ.If we choose the Frobenius norm ∥X ∥ F = Tr[X † X ], an explicit formula for χ + (t) exists [35].Since in our case χ(t) is Hermitian, Ref. [35] shows that the same expression is obtained by using the spectral norm ∥X ∥ ∞ = σ max (X ), where σ max (X ) is the maximum singular value of X .The result is that χ + (t) is the "positive part" of χ(t), obtained from χ(t) by putting to zero the negative eigenvalues: This gives a fairly general efficient way to determine χ + (t), at least numerically: it is sufficient to compute a spectral decomposition.
In order to gain some understanding we will now make some observations about the spectral structure of χ(t) in (8).For notational convenience, here we will not write the time parameter and we will use collective indices i = (k, q) and j = (n, m), lexicographically ordered.Moreover we write Γ αβ,i to mean Γ αβ (ω kq , t).Then we have where Up to now |A α 〉 and |G α 〉 are general vectors.Since every Hermitian matrix can be written in the form (20), it is quite difficult to say something general about its spectrum.
A case that can be treated explicitly is when there is only one noise channel M = 1.Here we can drop the α, β indices and obtain Let us ignore the trivial cases in which A i ≡ 0 or Γ i ≡ 0, which would lead to χ = 0. Then it is easy to see that the vectors |A〉 and |G〉 are linearly independent, unless Γ i ≡ Γ ̸ = 0, in which case |G〉 = Γ |A〉.In the latter scenario χ = 2 Re Γ |A〉〈A|: this is a rank-one matrix with nonzero eigenvalue λ = 2 Re Γ ∥A∥ 2 and associated normalized eigenvector |A〉 /∥A∥, where ∥A∥ 2 = i |A i | 2 .This case is not so interesting because, at least when t → ∞, λ ≥ 0 by Bochner's theorem and no regularization is needed.However note that many common models based on qubits and harmonic oscillators resonantly coupled with a bosonic bath fall exactly in the case mentioned above (see App. B for details).Suppose instead that |A〉 and |G〉 are independent.Then χ is a rank-two matrix and the eigenvectors associated with nonzero eigenvalues are of the form |v〉 = a |G〉 + b |A〉.Writing χ |v〉 = λ |v〉 and equating coefficients we find two solutions: Notice that by the Cauchy-Schwarz inequality λ + ≥ 0 and λ − ≤ 0, and we confirm that χ is not positive semidefinite in general [10], even in the time-dependent case.
It is instructive to rewrite these expressions in terms of physical quantities, which are encoded in Γ .Given a vector x ∈ N 2 let us define a notation that treats |A i | 2 as a probability distribution.Then where we defined for convenience the quantity and Var(x) := 〈x 2 〉 − 〈x〉 2 is the variance of x ∈ N 2 .A similar result was obtained in Refs.[15,25,34].In particular, in Ref. [25] the authors parametrically splitted the Redfield equation in a "positive" and a "negative" contribution, minimizing the latter with an optimized choice of the parameters.For a single noise channel this is equivalent to our formulation, since if we write χ = χ + + χ − , where χ − := (χ − χ † χ)/2, then ∥χ − χ + ∥ = ∥χ − ∥ and we know that χ + minimizes ∥χ − P∥ for positive semidefinite P. Our approach generalizes this view, since it provides a well-defined procedure to regularize the Kossakowski matrix with an arbitrary number of (even correlated) noise channels.
Notice that the bigger the variances in Eq. ( 25) the bigger the magnitude of the negative eigenvalue λ − , hence the regularization is expected to cause minimum disturbance when the environment correlation function is quite flat over the set of Bohr frequencies of the system.This is consistent with a Markovian dynamics requirement.In fact, the magnitude of the negative eigenvalue of the Kossakowski matrix has been used before to quantify non-Markovianity [36], and is related to the well-known non-Markovianity measure introduced by Rivas, Huelga, and Plenio in Ref. [37].
To conclude, let us provide the expression for the regularized Kossakowski matrix in the single noise channel scenario [cf.Eq. ( 21)].It is given by χ + = λ + |+〉〈+|, where |+〉 is the normalized eigenvector associated with λ + .For simplicity, let us indicate here λ ≡ λ + .Given the shape of a/b in (22), consider the eigenvector A calculation shows that The normalized eigenvector is therefore and then Using the expressions given above, one finds after some algebra that the components of χ + are 4 Example: open three-level system Now we want to compare different regularization procedures for a system that can be solved exactly, in order to try to assess which performs better.We will study a minimal model that is sufficiently complex to show the feasibility of our approach.A well-known example of an exactly solvable open quantum system is the spontaneous decay of a qubit into the field vacuum [6].As it will be clear later, this setting is not complex enough for our purposes since the Kossakowski matrix turns out to be rank-one positive in the limit t → ∞.A less known fact is that some kinds of open three-level systems can also be solved exactly [38], so we take that route.Consider the V-system depicted in Fig. 1, with free Hamiltonian (we assume the ground state |0〉 to be at zero energy).The environment is an infinite collection of bosonic modes with Hamiltonian H E = p ε p b † p b p and placed in the respective vacuum |Ω〉.The interaction is a linear coupling that causes transitions between the levels |0〉 and |1〉 and between the levels |0〉 and |2〉: where the coupling constants g 1,p and g 2,p are assumed for simplicity to be real numbers.
For the purpose of writing the Redfield equation, a simple calculation shows that the only relevant coupling operators are The others will make c αβ = 0 and hence do not appear in the final master equation.Instead, for these we have c αβ (τ) = p g α,p g β,p e −iε p τ ̸ = 0.For example, let us assume where µ, ω 0 > 0 and γ αβ = γ α γ β with γ 1 , γ 2 > 0. This exponential shape comes from a Lorentzian bath assumption with This is the choice that was made in Ref. [38] and we follow it here to provide a direct comparison between the exact dynamics and the various master equations.In order to make the paper self-contained we provide in App.C the derivation of the exact solution that is used in the following numerical calculations [cf.Eq. (60) and Eq. ( 69)].
In the simplifying case γ 1 = γ 2 = γ, we can consider γ as an estimate of the inverse evolution time of the system: γ ∼ 1/τ S .Moreover, we can take µ as an estimate of the inverse decay time of environment's correlations: µ ∼ 1/τ E .As a consequence, the Markovian approximation consists in assuming γ ≪ µ.

Numerical comparison
The structure of the Redfield equation [cf.Eq. ( 8)] is determined by the presence of the factor A β,kq A * α,nm = 〈k|0〉 〈β|q〉 〈0|n〉 〈m|α〉 .(35) This means that the only nonzero entries of χ(t) occur for k = n = 0 and q = β, m = α.With a quick calculation one realizes that where and the Lamb shift is We can also conveniently rewrite where ρ βα (t) = 〈β|ρ(t)|α〉 and where we dropped the time dependence for notational convenience.This is a linear system of first-order differential equations that can be efficiently solved by a numerical routine: here we adopted a Runge-Kutta method (RK45) provided by the Python library SciPy [39,40].Notice that the Kossakowski matrix for this setting (which is 9 × 9) is filled with zeros except for a 2 × 2 block with entries d αβ .Therefore it is clear that regularizing χ is equivalent to regularizing d.By choosing a two-level system instead of a three-level one the nonzero block would have consisted of a single entry: this scenario would be trivial since positivity is then guaranteed by Bochner's theorem, at least in the time-independent case.See App.B for clarifications on this point.Now we present a comparison between the exact solution provided by Ref. [38] [reported here in App.C, cf.Eq. (60) and Eq. ( 69)] and what we obtain by numerically solving the system in (41) for various choices of regularization of the matrix d.In Fig. 2 we report the results for the evolution starting from the pure initial state |ψ 0 〉 = (|1〉 + |2〉)/ 2, and choosing as parameters ω 1 = 1, ω 2 = 2, ω 0 = 1.5, γ 1 = γ 2 = 0.3, and µ = 2.At this level all equations behave more or less similarly.Except for the fact that the secular-approximated one globally provides the worst results, it is hard to tell which of the others performs better.The situation is similar for other choices of parameters.(14)] is obtained by finding the smallest coarse-graining time that guarantees positivity of the Kossakowski matrix.The label "ULE" refers to the "universal Lindblad equation" described in Ref. [22] [also, cf. ( 16)]. Here

Choi operator technique
If we want to assess more carefully the quality of a regularization procedure we should find a way to compare the results with the exact one in a way that is independent from the initial state.In order to do that, let us step back to the dynamical semigroup picture of Sec. 2. What we actually want is to compare the semigroup {Φ t,s } generated by our master equation with the semigroup {Φ (e) t,s } generated by the exact dynamics.Here we propose a simple approach to compare them pointwise, i.e., at fixed time t.More global comparisons should be possible but are out of the scope of the present paper and are left to future work.
Given a map Φ t,s : L(H S ) → L(H S ) we can construct the Choi operator [27, 28] A well-known fact is that Φ t,s is completely positive if and only if J (Φ t,s ) ≥ 0. However, we are mostly interested in the fact that there exists a bijection Φ t,s ↔ J (Φ t,s ), which is the Choi-Jamiołkowski isomorphism.The usefulness of this observation is twofold.First of all, on L(H S ⊗ H S ) we have well-established metrics that we can use, such as the Frobenius norm.
Secondly, if we compute we have a (pointwise) measure of the difference between the two dynamics that does not depend on the initial state.In Fig. 3 we report examples for δ(t) calculated with the Frobenius norm for various master equations.Here we fix ω 1 = 1, ω 2 = 2, ω 0 = 1.5, γ 1 = γ 2 = 0.05, and we present plots for several values of the spectral width µ.We find that the performance of our approach with respect to the others depends on the ratio between µ and the frequency range of the system ω R = max {ω 1 , ω 2 }.For values of µ sufficiently higher than ω R our procedure has little effect on the already good accuracy of the Redfield equation, as expected from the fact that the Kossakowski matrix is essentially positive in a deep Markovian regime (see discussion in Sec.3.2).For smaller values of µ (while still being greater than ω R ) we can instead observe how our version of the regularized Redfield equation approximates well the exact dynamics in a more consistent way with respect to the other master equations, especially at short times.This was not clear with a direct comparison at the level of the density matrix, and it is a consequence of our ability to retain time dependence in the Kossakowski matrix.
However, a change in the trend can be observed when lowering the value of µ below ω R , where our regularization provides results worse than Redfield itself.In this essentially non-Markovian regime the truncation of the negative part of the Kossakowski matrix has a too drastic effect on the dynamics.This situation is of course out of reach for the approach of the present paper, since it enforces CP-divisible dynamics, and other methods should be used.

Conclusions
In this work we looked at the problem of finding a LGKS-like equation from the microscopic dynamics as a regularization process of the Kossakowski matrix in the Redfield equation.With this picture in mind, we proposed to replace such a matrix with its closest positive semidefinite one, thus providing the CP-divisible dynamics that is closest to the Redfield one.We also used the Choi-Jamiołkowski isomorphism to envision a pointwise measure of the distance between two dynamical processes, and we applied it to the problem of assessing which master equation better approximates the exact dynamics of a simple open three-level system.We found our proposal to lead to the overall best results in this regard, provided one works in a Markovian regime where the spectral width of the environment is greater than the frequency range of the system.Notably, our approach is tailored to retain the time dependence of the Kossakowski matrix, allowing it to be accurate at short times.Unfortunately, at this level the approach is mainly numerical and it is still an open problem to understand what are the implications of the proposed manipulation on the thermodynamics of the system and the steady-state manifold structure.Note that at the Redfield level these kinds of characterizations already present some subtleties: for a system in contact with a finite-temperature bath the steady state is not the Gibbs state of the system and corrections should be included by considering a mean force Gibbs state [41].Such calculations depend on the shape of the Kossakowski matrix, which we are modifying in an analytically unpredictable way (at least in the general case), thus making an immediate transition to the regularized scenario difficult.
A possible future improvement would be to envision an alternative regularization scheme that is able to retain the non-Markovian features of the Redfield equation.Moreover, it would be desirable to have a meaningful measure of the distance between two dynamical processes which goes beyond the pointwise approach followed here: this is an interesting problem on its own and can lead to other general applications.
Another question that needs to be addressed is to what extent our conclusions can be applied to infinite-dimensional systems, where it is trickier to apply the LGKS theorem and where we expect the choice of the involved norms to matter more.
where we introduced the quantity Γ αβ defined in Eq. (10) of the main text.The exponential factor in (45) can be eliminated by going back to the Schrödinger picture: where we defined the matrix Let us focus on the term inside the round brackets in Eq. ( 46).If we expand the commutator and write explicitly the "H.c." part we get where we dropped the time dependence of ρ(t) and K kq,nm (t) for ease of notation.It is convenient to treat the first and third term of the right-hand side together.Replacing them in the sum in Eq. ( 46) we have k,q,n,m The second and fourth term of Eq. ( 48) can be treated similarly: Equation ( 7) is obtained by plugging Eqs. ( 49) and (50) in Eq. ( 46) and introducing the matrices η kq,nm , χ kq,nm defined in Sec. 2 of the main text.

B Kossakowski matrix for a qubit and a harmonic oscillator
In this section we compute the Kossakowski matrix of two common scenarios, in which a bosonic bath is coupled with either a qubit or a harmonic oscillator.Although simple, these examples share an interesting peculiarity: the Kossakowski matrix of the corresponding timeindependent Redfield equation is positive semidefinite, so that no regularization is needed to ensure the positivity of the dynamics.
In the first model a qubit is coupled to a bath of harmonic oscillators with a rotatingwave interaction Hamiltonian.Denoting the two energy levels of the qubit as |0〉 , |1〉 we have

while the interaction Hamiltonian writes
where b p is the creation operator relative to the environmental mode p.Using the notation of Eq. ( 4) we have A 1 = |0〉〈1|, and B 1 = p g 1,p b † p , while the Hermitian conjugates of A 1 , B 1 do not contribute to the dynamical equation since the relative correlation function 〈 b † p (τ)b p 〉 vanishes.The coordinates of A 1 in the decomposition (6) are simply given by A 1,kq = δ k,0 δ q,1 so that Then the Kossakowski matrix reads where we used ω 01 = ω 1 .The matrix in Eq. ( 53) is diagonal, with the single non-zero element being χ 01,01 (t).After the second Markov approximation is applied, we are left with χ 01,01 = 2 Re Γ 11 (ω 1 ) that is ensured to be positive as a consequence of Bochner's theorem.The case of the harmonic oscillator can be treated similarly.While the bath Hamiltonian is the same as the preceding example, we have H S = ω S a † a and where a, a † are creation and annihilation operators of the system.We repeat the calculations done for the qubit, but considering A 1 = a and obtaining We perform the sum on l 1 , l 2 and compute the associated Kossakowski matrix, that writes Looking for a redefinition of the indices (k, q) = i and (n, m) = j as in Sec.3.2, we notice that n, q are forced to be equal to m+1, k +1 respectively.The only ordered couples with a nonzero contribution to the r.h.s. of (56) are of the form (k, q) = (i, i + 1), so that we can adopt the simple mapping (k, q) = (i, i + 1) → j and (m, n) = ( j, j + 1) → i.

C Exact solution of the open three-level system
In this appendix we provide the exact solution of the open three-level system described in Sec. 4, which is originally described in Ref. [38].Suppose that the initial state of the universe is the following pure state: Since the total number of excitations is conserved, the state at time t must be of the form where α ∈ {1, 2}.From Eq. (61a) we immediately conclude that a 0 (t) = a 0 (0) for all t ≥ 0. Remembering that d p (0) = 0, Eq. (61c) can be formally integrated as where Q(s) = s 3 + h 1 s 2 + h 2 s + h 3 is a polynomial in s with coefficients Assuming that Q(s) has three non-degenerate roots r 1 , r 2 , r 3 we can apply the following Lagrange partial fraction decomposition: where Q ′ (s) = 3s 2 + 2h 1 s + h 2 is the derivative of Q(s).The result is â1 (s) = The inverse Laplace transform of these expressions leads to the desired solution: together with a 0 (t) = a 0 (0).

Figure 1 :
Figure 1: Open three-level V-system described in Sec. 4. The arrows indicate transitions induced by an external bosonic environment in its vacuum state.