Convergence guarantees for discrete mode approximations to non-Markovian quantum baths

Non-Markovian effects are important in modeling the behavior of open quantum systems arising in solid-state physics, quantum optics as well as in study of biological and chemical systems. The non-Markovian environment is often approximated by discrete bosonic modes, thus mapping it to a Lindbladian or Hamiltonian simulation problem. While systematic constructions of such modes have been previously proposed, the resulting approximation lacks rigorous and general convergence guarantees. In this letter, we show that under some physically motivated assumptions on the system-environment interaction, the finite-time dynamics of the non-Markovian open quantum system computed with a sufficiently large number of modes is guaranteed to converge to the true result. Furthermore, we show that this approximation error typically falls off polynomially with the number of modes. Our results lend rigor to classical and quantum algorithms for approximating non-Markovian dynamics.

Quantum systems invariably interact with their environment, and any simulation technique used to model their behaviour needs to capture this interaction. Traditionally, such interactions are analyzed within the Markovian approximation, wherein the system dynamics is described by the Lindbladian master equation [1]. However, a number of quantum systems arising in solid-state physics [2][3][4][5], quantum optics [6][7][8][9][10] as well as quantum biology and chemistry [11][12][13][14] cannot be modeled accurately within the Markovian approximation and the non-Markovian nature of the environment needs to be explicitly taken into account.
Simulating non-Markovian open quantum systems is difficult since it is usually not possible to formulate a dynamical equation for the system state from a given physical model of the system-environment interaction. While it is generically expected that non-Markovian open quantum systems satisfy a master equation of the Nakajima-Zwanzig form [13,15], or the time-convolutionless form [16][17][18], it is usually hard to obtain an explicit form of such a master equation except when the system is only weakly coupled to its environment [19][20][21][22]. Even though significant progress has been made in utilizing influence functionals or their variants for describing and simulating non-Markovian dynamics [23][24][25][26][27][28][29][30][31], the worst-case classical complexity of these approaches increases exponentially with time. An alternative approach is to identify and track a set of discrete modes that approximate the environment [32][33][34]. This maps the simulation of the non-Markovian open quantum system to a larger Hamiltonian or Lindbladian simulation problem, which can be solved using standard classical [35][36][37][38][39] or quantum algorithms [40][41][42][43][44][45].
For gaussian bosonic environments [46,47], there are two prominent approaches to identifying these discrete modes. First is to use the Lorentzian pseudomode theory [48][49][50][51][52][53], wherein the spectral density function of the non-Markovian environment is approximated by a finite sum of Lorentzians, each of which corresponds to an individual bosonic mode coupled to Markovian reservoir. The second method is to use * rahul.trivedi@mpq.mpg.de star-to-chain transformation [33,34], which uses the Lanczos iteration to identify a 1D chain of discrete bosonic modes with nearest neighbour couplings that approximate the environment and map the problem of computing non-Markovian quantum dynamics to a Hamiltonian simulation problem. While a systematic construction of these approximations have been given, their convergence properties are less well understood. For instance, classes of models where the Lorentzian pseudomode description is exact are known [48,49], but it is unknown if it can efficiently approximate a non-Markovian environment. More attention has been paid towards investigating the convergence properties of the star-to-chain transformation -convergence guarantees have been provided for models with a finite Lieb-Robinson velocity [54,55], or for bounded memory kernels [56] that can possibly have an infinite Lieb-Robinson velocity [57]. However, these analyses do not extend to models with distributional (and hence unbounded) memory kernels, such as those commonly encountered in non-Markovian quantum optical systems [6,58].
In this letter, we provide general and rigorous convergence guarantees for discrete-mode approximations of non-Markovian gaussian bosonic environments. We show that for a wide and physically-motivated class of non-Markovian models, both the Lorentzian pseudomode approximation and the star-to-chain transformation is guaranteed to converge and the approximation error falls off polynomially with the number of pseudomodes. Our results not only provide the first set of rigorous convergence guarantees for the pseudomode approximation, but also extend the convergence guarantees for the starto-chain transformation to non-Markovian systems described by a distributional memory kernel, such as those commonly encountered in quantum optics [6,58].
In order to perform this convergence study, there are several theoretical challenges that our work resolves. One of the key issues with analyzing models of non-Markovian environments is that the environment can support arbitrarily large energies and the quantum system often has non-vanishing couplings with high energy environment modes [6,58]. We identify a set of physically motivated sufficient mathematical conditions on the system-environment dynamics that allow for rigorously neglecting the high energy modes in the environment while studying finite-time dynamics. While this straightforwardly follows for problems where the system-environment couplings vanish at high energies, our analysis also includes distributional memory kernels. Combining this with an analysis of the Lorentzian pseudo-mode approximation and the star-to-chain transformation within a truncated environmentenergy window, we provide convergence guarantees as well estimates of the rate of convergence for both of these methods.
We consider an open quantum system model, with a ddimensional quantum system with hilbert space H S = C d (referred to as a 'local system') interacting with a gaussian environment whose Hilbert space, H E is assumed to be a fock space over L 2 (R). We denote by a ω the annihilation operator for this Fock space and consider Hamiltonians of the form where L : H S → H S is the operator describing the coupling of the system with the environment, and v is the frequency dependent coupling function between the environment and the system. We point out that the Hamiltonian in Eq. 1 is only a provably (essentially) self-adjoint operator if v ∈ L 2 (R) † . However, several problems in quantum optics (e.g. systems with point coupling) are described by v ∈ C ∞ b (R) ∩ S ′ (R) † † that are tempered distributions and still result in well defined local system dynamics. Throughout this paper, we will be interested in approximating the dynamics of the reduced state of the local system -for simplicity, we restrict ourselves to the case where the environment is initially in a vacuum state, extensions of the main results of this paper to initially excited environment states is provided in the supplement [59].
First we outline two physically motivated assumptions that we make on the model under consideration that are sufficient conditions to allow us to neglect large environment frequencies. The first assumption makes mathematically precise the expectation that the effective time-domain kernel, K v corresponding to the coupling function is approximable by its restriction to a finite frequency window. Since K can in general be a distribution, we first introduce a distributional quasi-norm to quantify this approximation error.
We note that our definition of this quasi-norm involves two applications of the kernel on a test function, which † A map v : is the space of tempered distributions that correspond to multiplication with bounded and smooth functions. itself is a function of two time indices -this choice of the quasi-norm is natural since we are interested in quantifying the back-action of particles emitted into the environment on the local system dynamics, with there being two time indices needed to describe the reduced state of each particle in the environment.
By ensuring that the kernel is distributionally approximable within a finite frequency window, assumption 1 ensures that the model described by Eq. 1 does not suffer from an ultraviolet divergence arising due to the environment being able to support arbitrarily high frequencies [60]. A number of models commonly considered in practice do satisfy this assumption satisfy assumption 1 (proposition 1 in the supplement [59]). Physically important examples of such environments include environments with a Lorentzian coupling function which is typically used to model an atomic system interacting with an optical cavity [61].
3. Environments modelling retardation effects, described by coupling functions of the form v(ω) = M i=1 v i e iωτi are commonly used to model retardation effects [62]. These environments also satisfy the conditions of assumption 1 (proposition 2 in the supplement [59]).
We point out that while the kernel might be approximable within a finite frequency window, due to the infinitedimensional nature of the environment's Hilbert space, it does not immediately follow that the reduced state of the local system is also similarly approximable. To proceed further, we introduce a second assumption that ensures physically meaningful joint system-environment dynamics. To state this assumption mathematically, we introduce the N −point Green's function of the localized system.
where U (τ 1 , τ 2 ) is the propagator corresponding to the Hamiltonian in Eq. 1 and L(τ ) = U (0, τ )LU (τ, 0) is the operator L in the Heisenberg picture. § For a Hilbert space H, L(H) is the set of operators T : H → H which are bounded. For T ∈ L(H S ), · is the usual operator norm defined by T = sup |ψ ∈H S \{0} T |ψ / |ψ . The physical significance of the N −point Green's function is that it determines the projection of the environment state on the N −particle subspace, as is made explicit in the following lemma (proved in the supplement [59]).
Lemma 1 Let |ψ(t) = U (t, 0) |σ ⊗ |vac , where |σ is a system state, then and ∀ω ∈ R N , Furthermore, as is made precise in lemma (proved in supplement [59]), the N −point Green's function is bounded, and consequently the norm of the N −particle projection of the environment state is also bounded.
Our second assumption can be interpreted as a bound on the rate at which the local system can emit or absorb an excitation from the environment. Any physically reasonable model of the environment, despite its non-negligible couplings with high frequency environment modes, is expected to satisfy this assumption.
This assumption can be proved for two cases: 1. Markovian environments, i.e. environments with a frequency independent coupling constant (v(ω) = v 0 ). In this case, an application of the quantum regression theorem can be used to show that assumption 2 is satisfied for such environments (proposition 3 in the supplement [59]).

2.
Environments with a square integrable coupling constant, assumption 2 can again be rigorously proven (proposition 4 in the supplement [59]).
With these two assumptions, we can now prove the convergence of the pseudomode theory [48] and star-to-chain transformation [33] for simulating non-Markovian quantum systems. Our first result, proved in the supplement [59], rigorously shows that in a finite amount of time, the localized system cannot excite arbitrarily high frequencies in the environment, and the environment can thus be approximated within a finite energy window. More precisely, is a coupling function such that assumptions 1 and 2 are satisfied. Denoting by ρ(t) the reduced density matrix of the local system at time t when an initial state |σ ⊗ |vac is evolved under the Hamiltonian in Eq. 1 and by ρ ωc (t) the reduced density matrix of the local system at time t when the same initial state is evolved under the Hamiltonian then where ε(ω c , t) is the cutoff error given by Here V (ω c , t) is defined in assumption 1 and where γ(t) is introduced in assumption 2.
The dependence of the cutoff error on ω c is a sum of two terms -one that falls off as O(1/ √ ω c ), which can be interpreted as the consequence of introducing a rectangular frequency window on the emitted photon wave-packet, and a second term which depends on the error introduced in approximating the time-domain kernel K v within a frequency window (assumption 1). The magnitude of this error, as given by the functions f 1 , f 2 , depends on the strength of the coupling between the system and environment (∝ v ∞ , L ), as well as on the rate at which particles are exchanged with the environment (γ(t) introduced in assumption 2). The result of theorem 1 is thus in-line with our intuition that a larger cutoff frequency is needed for systems which strongly couple to the environment, and absorb or emit particles very rapidly from the environment. Furthermore, the error grows exponentially with time t -this is a consequence of the fact that the local system can in principle emit an arbitrarily large number of particles into the environment while only being constrained by the bound in lemma 2. In practice, we expect the errors grow only polynomially with t, with the polynomial to depend on the maximum number of particles that the local system can emit into the en-vironment.
Next, we consider the pseudo-mode method [48], which approximates the non-Markovian dynamics of the local system by the Markovian dynamics of a larger system.
. . M }} corresponds to a coupling functionv which satisfies [48] |v(ω) In practice, to obtain a pseudomode description that approximates a given coupling function v, |v| 2 is approximated by a sum of lorentzians within a sufficiently large but finite frequency window, with each Lorentzian corresponding to an independent pseudomode. Our next result, proved in supplement [59], shows that this procedure is guaranteed to converge. The convergence rate of the pseudomode approximation will, in general, depend on the details of the coupling function v -in particular, on the growth of its derivative v ′ (ω) with ω, as well as on the falloff of the cutoff error ε(ω c , t), introduced in theorem 1, with the cutoff frequency. For typical coupling functions, |v ′ (ω)| = O(poly(ω)) and ε(ω c , t) = O(exp(O(t))poly(ω −1 c )). Under these assumptions, we show that the error incurred in the pseudomode approximation decreases polynomially with the number of pseudomodes.
is a coupling function such that assumptions 1 and 2 are satisfied and let ρ(t) be the reduced state of the local system after evolving an initial state |σ ⊗ |vac using the Hamiltonian in Eq. 1. Then, there exists a pseudomode description of the environment (definition 3) with M bosonic modes which provides an approximationρ(t) to the reduced state of the local system such that ρ(t) −ρ(t) tr → 0 as M → ∞. Furthermore, if |v ′ (ω)| = O(poly(ω)) and the cutoff error ε(ω c , t) = O(exp(O(t))poly(ω −1 c )), then there exists a pseudomode description of the non-Markovian system with M bosonic modes such that the trace-norm error in approximating the reduced local system state at time t scales as O(exp(O(t))poly(M −1 )).
Finally, we consider the star-to-chain transformation, which maps the non-Markovian environment to a 1D bosonic tightbinding model with the local system effectively coupled to the first bosonic mode.
with the environment through the operator L ∈ L(C d ) with the initial system-environment state being |σ ⊗ |vac , its reduced state at time t is given by ρ(t) = Tr aux (Û (t, 0) |σ σ|⊗ (|0 0|) ⊗MÛ (0, t)), whereÛ (τ 1 , τ 2 ) is the propagator corresponding to the Hamiltonian A chain transformation of the environment can be explicitly constructed by using the Lanczos iteration -this proceeds by first introducing a frequency cutoff ω c , and then starting from the mode a 0 ∝ ωc −ωc v(ω)a ω dω, applying the Lanczos iteration with respect to the environment Hamiltonian. This yields the parameters {(ω i , g i ) : i ∈ {0, 1, 2 . . . M − 1}} for the chain description of the environment. By exploiting the previously established connection between the star-to-chain transformation and orthogonal polynomials [33,63], we show (supplement [59]) that, under assumptions 1 and 2, the star-tochain transformation will converge for sufficiently large number of modes and that the convergence rate is polynomial in the number of modes.
is a coupling function such that assumptions 1 and 2 are satisfied and let ρ(t) be the reduced state of the local system after evolving an initial state |σ ⊗ |vac using the Hamiltonian in Eq. 1

. Then, there exists a chain description of the environment with M modes (definition 4) that provides an approximationρ(t) to the reduced state of the local system
, then the trace-norm error in approximating the reduced local system state at time t scales as O(exp(O(t))poly(M −1 )).
In conclusion, our work provides a rigorous analysis of Markovian dilations to non-Markovian open quantum systems. We show that the finite-time dynamics of a wide class of non-Markovian quantum systems can always be well approximated by a larger Markovian system, and we also provide theoretical scalings of how much larger the effective Markovian system is. Several questions of interest to open quantum system theory are left open in this work. One of the questions that we leave unanswered is to study the class of coupling functions, with possibly distributional kernels, for which the resulting system-environment dynamics is well defined. A rigorous study of this problem would be relevant to advancing the mathematical understanding of non-Markovian open quantum system models, and could lead to a general proof of assumption 2. Another direction to pursue would be to improve the exponential dependence of the error estimates on time t, or to identify settings in which these error estimates are tight.

Function spaces and analysis:
For v ∈ R N and u ∈ R M , we will denote by {v, u} a vector in R N +M whose first N entries are elements of v and the next M entries are elements of u.
All the integrals in this paper are Lebesgue integrals and with respect to the Lebesgue measure over R N . For a measure space (Ω, F, µ), given a measurable map f : Ω → R, we define its essential supremum via esssup x∈Ω f (x) = inf{a ∈ R | µ(f −1 ([a, ∞])) = 0}. The essential supremum can be interpreted as the supremum of a function obtained on ignoring subsets of its domain with zero measure. A propositional function P : Ω → {True, False} is said to be True almost everywhere in Ω if the µ {x | P(x) is False} = 0. This will be denoted as ∀x ∈ a.e. Ω : P(x) is True.
We will denote by L 2 (R D ) the set of complex-valued square integrable functions over R D . We will denote by B(R D ) the set of bounded functions over R D . For f ∈ B(R D ), we denote by f ∞ the maximum magnitude of the function i.e.
For a closed interval I ⊂ R, we will denote by C k (I) the set of functions f : I → C which are k−continuously differentiable. In particular, C 0 (I) is the set of continuous complex-valued functions with domain I and C ∞ (I) is the space of smooth complex-valued functions with domain I. We will also denote by C k b (I) the space of bounded k−times continuously differentiable functions -on compact intervals I, this will coincide with C k (I). We will denote the space of tempered distributions over R by is the space of tempered distributions that correspond to multiplication with bounded and smooth functions. Finally, for any α > 0, rect α : R → R to be denote the rectangular windowing function defined by For a < b, we will denote by AC([a, b]) the space of all functions f : [a, b] → C that are absolutely continuous. We recall An immediate consequence of this definition is that the function f has a derivative, which we denote by ∂ the space of symmetric, continuous functions f : [a, b] 2 → C which satisfy the following two conditions 1. f is absolutely continuous in either of its arguments i.e. ∀y ∈ [a, b], the map ϕ y : We note that the partial derivatives of f with respect to either and both of its arguments, A symmetric map K is also positive if The trace-norm will be used to quantify distance between two density operators. It is useful to note that given two states |ψ 1 , |ψ 2 ∈ H with unit norm, then the operator |ψ 1 ψ 1 | − |ψ 2 ψ 2 | is bounded and trace-class and its trace-norm satisfies This property, together with the contractivity of the trace-norm under partial trace, allows us to translate an error in between two states of a composite quantum system to an error between the density operators of a subsystem of the composite quantum system.

Error analysis:
Given two functions f, g :

II. N -POINT GREEN'S FUNCTION AND THE ENVIRONMENT STATE
In this section, we relate the N −point Green's function introduced in the main text to the joint state of the local system and environment (lemma 1) and also provide a short proof of lemma 2 from the main text. For generality, in this section we consider the case where the environment is in an initially excited state |E ∈ H E -the results of the main text correspond to the setting where |E ≡ |vac , while in section V, we consider an initially excited environment state.
is the propagator for the joint system-environment dynamics and L(τ ) = U (0, τ )LU (τ, 0) is the operator L in the Heisenberg picture.
where |σ ∈ H S is an initial system state and |E ∈ H E is the initial environment state, then Proof : The state |ψ(t) admits a representation of the form where . . d}} is a basis for the d−dimensional local system and where a ω (t) = U (0, t)a ω U (t, 0) and the time-ordering operator can be trivially introduced since all the operators are evaluated at the same time. From Heisenberg equations of motion, it immediately follows that from which it immediately follows that Substituting Eq. 5 into Eq. 4 and using d i=1 |e i e i | = I, we obtain that Proof : This follows from noting that Next, we note from Eq. 6 that This finishes the proof of the lemma. .

III. SETTINGS WHERE ASSUMPTIONS 1 AND 2 ARE PROVABLE
A. Assumption 1 For completeness, we restate assumption 1 from the main text, and then prove it for several classes of coupling functions encountered in practice.
Proposition 1 Assumption 1 is true for coupling functions v that are bounded and square integrable, v ∈ B(R) ∩ L 2 (R).
Proof : We consider constants v : R → C that are smooth, bounded and square integrable. We note that in this case, for a symmetric and positive kernel K : 2 , it follows that square-integrable coupling constants satisfy the conditions of assumption 1 with V (ω c , t) being given by We note that the case of a square integrable coupling constant contains the physically important case of a lorentzian coupling function, which often arises in studying the interaction of quantum systems with lossy resonators. Next, we show that assumption 1 is also true for environments whose coupling function is a sum of a discrete number of harmonic terms. Such coupling functions arise in models that consider a local quantum system with point couplings to a one-dimensional propagating field. We first provide a bound on the error incurred in approximating a delta function by a sinc function.
Proof : We begin by using the Dirchlet integral, ∞ 0 sin x/xdx = π/2, from which it follows that The second term is easily bounded To bound the first term, we split the integral at b(p) ∈ [0, a], where b(p), to be specified later, vanishes as p → ∞: Futhermore, using integration by parts, we obtain Using these estimates and noting that sup Since this bound holds for any b(p) ∈ [0, a], we choose b(p) = Θ( a/p) which proves the lemma .
Proof : This follows straightforwardly from integration by parts. sin px from which the lemma statement follows .
where ξ(a, b, y, p) → 0 as p → ∞ and δ(·) is the dirac-delta function, whose action on a continuous function f within an 2. If y < a or y > b, then it follows from lemma 2 that Here we have used lemma 1 and the fact that ∀y ∈ (a, b) .
Since ϕ y is absolutely continuous, it follows from lemma 3 that b a πδ( Furthermore, it also follows from lemma 3 and the fact that K was absolutely continuous in either of its arguments that Similarly, for x almost everywhere in [a, b], ϕ ′ y2 (x) exists and consequently from lemma 3 it follows that Using this along with Eq. 7 completes the proof .
V i e iωτi ∀ω ∈ R, then assumption 1 is satisfied for v.
Proof : We note that it is possible to express |v| 2 as a finite sum of harmonics: For any K ∈ AC sym ([0, t] 2 ), we obtain Applying lemma 4, we immediately obtain which proves that assumption 1 is satisfied . We point out that proposition 2 shows that assumption 1 is satisfied for both the Markovian model as well as models that model retardation effects.

B. Assumption 2
In this appendix, we consider two classes of open quantum systems wherein assumption 2 can be proven rigorously. The first is the case of a Markovian coupling function, and the second is the a non-Markovian coupling constant with a bounded square-integrable coupling function. Noting that To proceed further, we provide a bound on the second term in the above estimate. Denote by |φ = U (s, 0)L(s i )L(s i−1 ) . . . L(s 1 ) |σ ⊗ |vac , we can write |φ as a superposition of its N −particle components: where {|e i , i ∈ {1, 2 . . . D}} is a basis for the D−dimensional local system and where a ω (s) = U (0, s)a ω U (s, 0) and we have used s 0 ≥ s i ≥ s i−1 · · · ≥ s 1 to introduce the time-ordering operator. We then use the Heisenberg equations of motion for a ω (s) to express We thus obtain This immediately allows us to bound the norm of the M −particle component of |φ : Finally, noting that the restriction of the operator and consequently we obtain that ∀σ ∈ H S , s 0 ∈ (s i , s i+1 ), Since this bound holds for every i ∈ {0, 1 . . . N − 2}, it proves that assumption 2 holds for square integrable coupling functions .

IV. DETAILED PROOFS
A. Proof of theorem 1: Introducing a frequency cut-off on the environment We first show rigorously that, subject to the assumptions 1 and 2, a frequency cutoff can be introduced in the environment without significantly impacting the dynamics of the local system.

Definition 2 For a given frequency cutoff
Furthermore, P ωc , Q ωc : H S ⊗ H E → H S ⊗ H E such that P ωc = id ⊗ Π ωc and Q ωc = id − P ωc .
The operator P ωc thus projects a given system-environment state onto the space of states with no environment frequencies outside of |ω| ≤ ω c . It is straightforward to see that the Hamiltonian H ωc (t) defined in theorem 1 satisfy H ωc (t) = P ωc H(t)P ωc . We now study the error introduced by evolving an initial state under the Hamiltonian H ωc (t) instead of the Hamiltonian H(t). This is made explicit in the following lemma.
Lemma 5 Given an initial state |ψ 0 , and denoting by U (t 1 , t 2 ) and U ωc (t 1 , t 2 ) the propagators corresponding to H(t) and H ωc (t) respectively, the states |ψ(t) = U (t, 0) |ψ 0 and |ψ ωc (t) = U ωc (t, 0)P ωc |ψ 0 satisfy |ψ(t) − |ψ ωc (t) ≤ Q ωc |ψ 0 + Q ωc |ψ(t) + t 0 P ωc H(s)Q ωc |ψ(s) ds Proof : Note that |ψ(t) = P ωc |ψ(t) + Q ωc |ψ(t) . From the Schroedinger's equation for |ψ(t) , it then follows that which can be integrated to obtain P ωc |ψ(t) = U ωc (t, 0)Q ωc |ψ 0 + |ψ ωc (t) − i The two terms in the error estimate provided in this lemma can be respectively interpreted as the error introduced in quantum state of the system on ignoring the high-frequency components of the environment state, and the error introduced by not accounting for the coupling of the high-frequency components of the environment state to the low-frequency components by virtue of their interaction with the local system. In the following lemmas, we show that both of these errors are bounded by a function that goes to 0 as the frequency cutoff ω c is made increasingly larger.

Lemma 6
If assumptions 1 and 2 are satisfied, it follows that ∀t > 0, where γ(t) is introduced in assumption 2.
Proof : It follows from the definition of Q ωc that Using Eq. 6, we obtain Therefore, where, for N ≥ 1 and s ∈ [0, t] N −1 , H N (ω, s; t) is given by Using integration by parts, together with assumption 2, we immediately obtain that and therefore from which the lemma immediately follows .
Lemma 7 If assumptions 1 and 2 are satisfied, it follows that ∀t > 0, where γ(t) is introduced in assumption 2.
Proof : Using the definitions of P ωc and Q ωc , it follows that Therefore, To proceed further, we make use of assumption 1. For |σ . An application of the assumption 1 using this kernel yields ∀s ∈ [0, t] N −1 where V (ω c , t) is the function introduced in assumption 1. Using assumption 2, we immediately obtain that ∀s It thus follows that Using this bound, the lemma statement follows.
is a coupling function such that assumptions 1 and 2 are satisfied. Denoting by ρ(t) the reduced density matrix of the local system at time t when an initial state |σ ⊗ |vac is evolved under the Hamiltonian in Eq. 1 and by ρ ωc (t) the reduced density matrix of the local system at time t when the same initial state is evolved under the Hamiltonian then where ε(ω c , t) is the cutoff error given by Here V (ω c , t) is defined in assumption 1 and where γ(t) is introduced in assumption 2.
Proof : Let |σ ∈ H S be an arbitrary local system state, and denote by |ψ(t) the joint state of the system and the environment obtained on evolving |σ ⊗ |vac under the Hamiltonian H(t). Similarly, we denote by |ψ ωc (t) the joint state of the system and environment obtained on evolving |σ ⊗ |vac under the Hamiltonian H ωc (t) = P ωc H(t)P ωc . We denote by ρ(t) and ρ ωc (t) the reduced state of the local system at time t i.e. ρ(t) = Tr HB [|ψ(t) ψ(t)|] and ρ ωc (t) = Tr HB [|ψ(t) ψ(t)|]. We can then use the contractivity of the partial trace, as well as the relationship between the trace-norm distance and the norm distance between two pure states (Eq. 3), we obtain An application of lemmas 5, V and V allows us to bound |ψ(t) − |ψ ωc (t) and obtain the result in theorem .
and R(0) = |σ σ| ⊗ (|0 0|) ⊗M . We first develop bounds on the error incurred on approximating one bounded and square integrable coupling constant with another one which are precisely stated in the following lemma.

Lemma 8 (Ref. [3]) An environment interacting with a local system with a coupling function v which satisfies
be coupling functions and let H 1 (t), H 2 (t) denote the Hamiltonians corresponding to these coupling functions as given by Eq. 1 in the main text. Denoting by |ψ i (t) as the joint state of the local system and environment when an initial state |σ ⊗ |vac , with |σ ∈ H S , is evolved under the Hamiltonian H i (t), i ∈ {1, 2}, then where , quantifies a distance between the two coupling functions.
Proof : For i ∈ {1, 2}, we denote by U i (τ 1 , τ 2 ) the propagator corresponding to the Hamiltonian H i . Furthermore, let |ψ i (t) = U i (t, 0) |σ ⊗ |vac , where |σ ∈ H S . It can be seen that An integral equation can easily be derived by differentiating this equation, followed by integrating it from 0 to t to obtain from which we obtain the bound We note from the definition of the Hamiltonians Therefore, Note that the operators ∆ ± are unbounded, and consequently the norms in the above equations cannot be estimated trivially by introducing the operator norms of ∆ ± . However, ∆ ± are bounded when restricted to only the N −particle subspaces of the environment's Hilbert space. In particular, for a N −particle state |ψ N Furthermore, the N −particle component of the state |ψ 1 (t) has a bounded norm, with Eq. 5 of the main text providing the explicit bound. It thus immediately follows that These estimates, along with Eq. 14, proves the statement of the lemma. . While this estimate allows for an assessment of the error incurred in approximating one coupling function with another, if only the state of the local system is of interest, the phase of the coupling function is irrelevant. This is easily seen from the fact that the transformation a ω → a ω e iθ(ω) and v(ω) → v(ω)e −iθ(ω) , for any θ : R → R, leaves both the Hamiltonian in Eq. 1 of the main text and the initial state (provided the environment is in the vacuum state) unchanged. Consequently, a bound can be provided on the error in the reduced density matrix of the local system that only depends on the magnitude of the coupling constant, as made explicit by the following lemma.
be coupling functions and let H 1 (t), H 2 (t) denote the Hamiltonians corresponding to these coupling functions as given by Eq. 1 of the main text. Denoting by ρ i (t) as the reduced state of the local system when an initial state |σ ⊗ |vac , with |σ ∈ H S , is evolved under the Hamiltonian H i (t), i ∈ {1, 2}, then where g(t) is the function defined in lemma 9.
Proof : Denote by θ i (ω) the phase of the coupling coefficient v i (ω) ∀ω ∈ R, i ∈ {1, 2}. Applying the transformation a ω → a ω e −iθi(ω) to the Hamiltonian H i , we obtaiñ For i ∈ {1, 2}, denote by |ψ i (t) the state obtained by evolving |σ ⊗ |vac under the HamiltonianH i (t). While these states, in their entirety, will be different from the ones obtained on evolving the system under the Hamiltonian H i (t), i ∈ {1, 2}, the reduced density matrix of the local system will be the same. Consequently, we can bound the error in the reduced density matrix of the local system by simply bounding the error between the states |ψ 1 (t) and |ψ 2 (t) , which can be done with an application of lemma 9: where g(t) and d(v 1 , v 2 ) are defined in lemma 9. Furthermore, we note that and thus we obtain Finally, using Eq. 3 together with the fact that partial trace is contractive, we obtain the bound provided in the lemma .
To use pseudomodes to simulate the dynamics of a non-Markovian quantum system, we would need to approximate the magnitude-square of the coupling function, |v(ω)| 2 , by a sum of Lorentzian. Since, as shown in theorem 1, the dynamics of local system is only impacted by the coupling function within a finite frequency window, a natural procedure to obtain the sum of Lorentzian decomposition would be to choose a frequency cutoff ω c > 0, an integer M indicating the number of Lorentzians and solve the following optimization problem: While in practice this optimization problem can be solved numerically, a theoretical analysis of the resulting approximation error as quantified by the solution of this problem is hard since it is a non-convex problem. However, as is made concrete in the following lemma, an upper bound on the approximation error can be provided by an explicit Lorentzian construction.
Proof : Our goal is to develop a Lorentzian approximation to a function Γ which is bounded and differentiable. We develop such an approximation in two steps: 1. We first use the fact that Lorentzians are mollifiers to construct an approximation Γ 1 : R → R to Γ within the frequency Note that as κ → 0, Γ 1 (ω) is expected to converge to Γ(ω) for ω ∈ (−ω c , ω c ). For a non-zero value of κ, an error will be incurred in approximating Γ by Γ 1 .
2. Next, we numerically discretize the integral in Eq. 15 to obtain the approximationΓ : R → R to Γ. Denoting by the number of Lorentzians (i.e. number of points used in discretizing the integral) by M , the discretization used in approximating the integral is δ = 2ω c /(M − 1) and thusΓ We now consider the error incurred in approximating Γ within a frequency window [−ω c , ω c ] byΓ. The relevant error of interest is the L 1 norm error between Γ · rect ωc andΓ, which can be upper bounded by Error between Γ and Γ1 Error between Γ1 andΓ + |ω|≥ωc
We can estimate each of these terms individually.
Error between Γ and Γ 1 : It follows from the definition of Γ 1 that The first integral can analytically integrated to obtain: To estimate the second integral, we make use of Taylor's theorem: We thus obtain that ωc −ωc We note that ξ, ξ d = O(xlog(1/x)) and therefore ωc −ωc Error between Γ 1 andΓ: Next, to obtainΓ from Γ 1 , we simply discretize the integral in the definition of Γ 1 -the Taylor's theorem can again be used to bound the resulting discretization error. In particular, Error due to tails: This error can be explicitly evaluated.
Noting that for x > 0, tan −1 (1/x) is convex upwards, it follows that for any θ > 0 Consequently, This proves the lemma .
Theorem 2 (Pseudomode approximation, restated) Suppose v ∈ C ∞ b (R) is a coupling function such that assumptions 1 and 2 are satisfied and let ρ(t) be the reduced state of the local system after evolving an initial state |σ ⊗ |vac using the Hamiltonian in Eq. 1. Then, there exists a pseudomode description of the environment (definition 3) with M bosonic modes which provides an approximationρ(t) to the reduced state of the local system such that ρ(t) −ρ(t) tr → 0 as M → ∞. Furthermore, if |v ′ (ω)| = O(poly(ω)) and the cutoff error ε(ω c , t) = O(exp(O(t))poly(ω −1 c )), then there exists a pseudomode description of the non-Markovian system with M bosonic modes such that the trace-norm error in approximating the reduced local system state at time t scales as O(exp(O(t))poly(M −1 )).
Proof : For a coupling function v ∈ C ∞ b (R), we denote by ρ(t) the true reduced state of the local system at time t i.e. the reduced state obtained when evolved as per the Hamiltonian in Eq. 1 of the main text. First we introduce a frequency cutoff ω c -we denote by ρ ωc (t) the reduced state of the local system obtained when evolved as per Hamiltonian in Eq. 1 of the main text but with coupling function v · rect ωc . Second, we approximate v · rect ωc by a sum of lorentzians by using the construction of lemma 11 to obtain a pseudomode approximation to the non-Markovian system (lemma 8). Letρ(t) be the reduced state of the local system at time t obtained on using the pseudomode approximation. We estimate the error between the true reduced state ρ(t) and the reduced state using the pseudomode approximationρ(t) via Lorentzian approximation error . We note that since by definition Γ d max (ω c ) is an increasing non-negative function, is an increasing function, Using these choices of κ and ω c , the error in approximating |v| 2 with a sum of lorentzian, as calculated in Lemma 11, is O(M −1/8 ) and thus goes to 0 as M → ∞. Furthermore, since ω c grows with M , it follows from theorem 1 that the cutoff error goes to 0 as M → ∞.
To show that, under the assumption of a polynomial growth in |v ′ (ω)| with ω and a polynomial fall off of the cutoff error with ω, the total approximation error decreases polynomially with the number of pseudomodes, we simply have to note that in the above proof, f −1 5 (M 1/5 ) will be poly(M ) and consequently the construction of theorem 3 yields a cutoff frequency ω c which grows polynomially with M . It is already shown in the proof of theorem 3 that the error in approximating the magnitude square of the coupling function with a sum of Lorentzian reduces polynomially with the number of pseudomodes, and with this choice of ω c the cutoff error also reduces polynomially with the number of pseudomodes .

ALGORITHM:
• Initialize p 0 (ω) = 1, β 0 = 0. We make two remarks about the Lanczos iteration described in the above algorithm. First, the cutoff frequency ω c is introduced into the algorithm so as to ensure that the integrals involved in computing the chain parameters are well defined. This is not necessary if the coupling functions decays fast enough with frequency so as to ensure that the integrals exist [5]. Second, the chain transformation generated by the Lanczos iteration with M → ∞ exactly reproduces the dynamics of the local system interacting with an environment restricted to [−ω c , ω c ]. However, only using a finite number of modes results in a truncation error -in the following lemma, we first analyze this truncation error and then combine it with the error incurred with the introduction of the frequency cutoff.
Lemma 12 Let ρ M→∞ (t) and ρ M (t) be the reduced state of the local system obtained on evolving an initial state |σ ⊗ |vac , with |σ ∈ H S , with respect toĤ M→∞ (t) andĤ M (t) (defined in definition 4) respectively, then and with the initial state being |σ ⊗ |vac . We note that the norm of the N −particle component of |ψ M→∞ (t) will satisfy the bound in lemma 2 of the main text since it corresponds to the quantum state obtained from a Hamiltonian of the form Eq. 1 of the main text with coupling function v · rect ωc , and translating this state into the interaction picture does not change these norms. Consequently, we obtain Noting that for an N −particle wave-packet |ψ N and for all s ∈ [0, t] Using this together with the bound in lemma 2 of the main text, we obtain that Finally, the lemma follows using the contractivity of the partial trace together with Eq. 3 .
An explicit dependence of the error bound in lemma 12 on the number of modes N can be obtained by utilizing the connection between orthogonal polynomials and the gauss quadrature method for approximating integrals with summations. This is provided in the following lemma.

Lemma 13
Let δB M (s) be as defined in lemma 12, then Proof : We start by noting that H E,M→∞ = ωc −ωc ωa † ω a ω dω and therefore Evaluating g 0 e isHE,M a 0 e −isHE,M is more involved since we need to diagonalize H E,M . It is convenient to introduce the polynomial π i , which is p i (defined in the algorithm 1) after normalization We note that from the orthogonality of the eigenmodes φ i , it follows that Finally, consider evaluating e isHE,M a 0 e −isHE,M . Noting that a 0 = Next we use the Gauss quadrature theorem [6] -we note that the polynomials p 0 , p 1 , p 2 . . . generated in algorithm 1 are the same polynomials that would be generated if a function of the form P (ω)|v(ω)| 2 for some P has to be integrated in the interval [−ω c , ω c ] using Gaussian quadrature. The M roots of p M correspond to the points within the interval [−ω c , ω c ] at which the integrand needs to be evaluated to compute the integral. It is well known that the gaussian quadrature with M such points exactly evaluates the integral if P is a polynomial of degree ≤ 2M − 1. More specifically, for a given M > 0, ∃w 0 , w 1 . . . w M−1 > 0 with M−1 i=0 w i = 1 such that for all polynomials P of degree ≤ 2M − 1, evolved as per Hamiltonian in Eq. 1 but with a frequency cutoff. Second, we use the Wilson's star-to-delta transformation with M modes to obtain an approximationρ(t) to ρ ωc (t). We estimate the error between the true reduced state ρ(t) and the reduced state using the pseudomode approximationρ(t) via ρ(t) −ρ(t) tr ≤ ρ(t) − ρ ωc (t) tr Cutoff error + ρ ωc (t) −ρ(t) tr Wilson star-to-delta error .
From theorem 1, it follows that the cutoff error decreases as ω c → ∞. To estimate ρ ωc (t) −ρ(t) tr , we use lemmas 12 and 13, which allows us to relate this error to the frequency cutoff ω c and the number of modes M . We note that choosing ω c to increase sublinearly with M , both the cutoff error ρ(t) − ρ ωc (t) tr and ρ ωc (t) −ρ(t) tr go to 0 as M → ∞, thus providing a convergence guarantee for Wilson star-to-delta transformation. Furthermore, with this choice of ω c and under the assumption that the cutoff error falls polynomially with ω c , it immediately follows that the approximation error decreases polynomially with the number of modes M .

V. THEOREM 1 FOR AN INITIALLY EXCITED ENVIRONMENT
In this section, we consider the setting where the environment is initially in an excited state. Some minor modifications to the formulations of assumptions 1 and 2 allows us to extend theorem 1 to an initially excited environment state. We consider an initially mixed environment state, described as an ensemble ρ E (0) := {p(E), |E ∈ H E }. We first assume that this ensemble is itself approximable within a finite frequency window, as is made precise below. where Π ωc is defined in definition 2.
Furthermore, a physically reasonable environment state should modify assumption 2 to bound the derivative of Green's function over the pure states in the ensemble describing the mixed environment state. Theorem 1 ′ Suppose v ∈ C ∞ b (R) is a coupling function and ρ E (0) is an initial environment state such that assumptions 1, 2 ′ and 3 are satisfied. Denoting by ρ S (t) the reduced density matrix of the local system at time t when an initial state ρ S (0)⊗ ρ E (0) is evolved under the Hamiltonian in Eq. 1 of the main text and by ρ S,ωc (t) the reduced density matrix of the local system at time t when the same initial state is evolved under the Hamiltonian then ρ(t) − ρ ωc (t) tr ≤ ε(ω c , t), where ε(ω c , t) is the cutoff error given by Here V (ω c , t) is defined in assumption 1 and where γ(t) is introduced in assumption 2 ′ .
Proof : Suppose that ρ E (0) = {p(|E ), |E ∈ H E }, is the initial environment state, and suppose ρ S (0) = d i=1 p i |σ i σ i | for some |σ 1 , |σ 2 . . . |σ d ∈ H S . The initial joint state of the system and the environment can then be expressed as p i dp(|E ) |σ i σ i | ⊗ |E E| .
We will first analyze the cutoff error obtained with an initial state |σ i ⊗ |E for some i ∈ {1, 2 . . . d} and |E ∼ ρ E (0). We note that for this problem, lemmas 5, and are applicable if assumptions 1, 2 ′ and 3 hold, and we obtain that holds for all i ∈ {1, 2 . . . d} and for almost all |E drawn from the ensemble ρ E (0). Thus, we immediately obtain that From Eqs. 22 and 23, the lemma statement follows. .