# Convergence of numerical approximations to non-Markovian bosonic gaussian
environments

Rahul Trivedi1,2<EMAIL_ADDRESS>Daniel Malz1,2 J. Ignacio Cirac1,2
1Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Str. 1, 85748 Garching,
Germany.
2Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4,
D-80799 Munich, Germany.

###### Abstract

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. A common approach to the analysis of such
systems is to approximate the non-Markovian environment by discrete bosonic
modes thus mapping it to a Lindbladian or Hamiltonian simulation problem.
While systematic constructions of such modes have been proposed in previous
works [_D. Tamascelli et al, PRL (2012), A. W. Chin et al, J. of Math. Phys
(2010)_], the resulting approximation lacks rigorous convergence guarantees.
In this paper, we initiate a rigorous study of the convergence properties of
these methods. 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, for most
physically interesting models of non-Markovian environments, the approximation
error falls off polynomially with the number of modes. Our results lend rigor
to numerical methods used for approximating non-Markovian quantum dynamics and
allow for a quantitative assessment of classical as well as quantum algorithms
in simulating non-Markovian quantum systems.

##  I Introduction

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 Breuer  _et al._ (2002). However, a number of quantum systems
arising in solid-state physics Finsterhölzl  _et al._ (2020); Chin  _et al._
(2011); Groeblacher  _et al._ (2015); de Vega  _et al._ (2008), quantum optics
Calajó  _et al._ (2019); Andersson  _et al._ (2019); González-Tudela  _et al._
(2019); Aref  _et al._ (2016); Leonforte  _et al._ (2021) as well as quantum
biology and chemistry Ishizaki and Tanimura (2005); Chin  _et al._ (2012);
Ivanov and Breuer (2015); Caycedo-Soler  _et al._ (2021) cannot be modelled
accurately within the Markovian approximation and the non-Markovian nature of
the environment needs to be explicitly taken into account.

One of the key difficulties with simulating non-Markovian open quantum systems
is that it is usually not possible to write a master equation that captures
the dynamics of this system from a given physical model of the system-
environment interactions. While it is generically expected that non-Markovian
open quantum systems satisfy a master equation of the Nakajima-Zwanzig form Xu
_et al._ (2018); Ivanov and Breuer (2015), or the time-convolutionless form
Smirne and Vacchini (2010); Pereverzev and Bittner (2006); Kidon  _et al._
(2015), it is usually hard to obtain an explicit form of such a master
equation except for when the system is only weakly coupled to its environment
in which case a perturbative treatment is possible Vacchini and Breuer (2010);
Schröder  _et al._ (2007); Timm (2011); Mukamel  _et al._ (1978). Another line
of work considers the time-evolution of non-Markovian open quantum systems as
the process of updating a tensor on the joint state-space of the system and
its past copies Albeverio  _et al._ (2007); Makarov and Makri (1994); Smith
and Caldeira (1987); Whalen  _et al._ (2017); Grimsmo (2015). Several
techniques from tensor networks can be imported to efficiently performing this
update, which allows for a numerically exact simulation of the non-Markovian
effects in the system Jørgensen and Pollock (2019); Del Pino  _et al._ (2018);
Strathearn  _et al._ (2018); Rosenbach  _et al._ (2016).

An alternative approach to simulating the dynamics of the non-Markovian open
quantum systems is to keep track of the state of the environment Huh  _et al._
(2014); Chin  _et al._ (2010); Woods  _et al._ (2014); Fang (2019). Since
keeping track of the entire environment state is computationally inefficient
in most cases, these methods start by identifying the most important modes in
the environment and reducing the simulation of the non-Markovian open quantum
system to a Hamiltonian simulation or Lindbladian simulation problem. The
resulting problem can then be efficiently solved using standard classical
Vidal (2004); Daley  _et al._ (2004); Banuls  _et al._ (2009); Daley  _et al._
(2009); Verstraete  _et al._ (2004) or quantum algorithms Berry  _et al._
(2014); Low and Chuang (2019); Berry  _et al._ (2015); Zanardi  _et al._
(2016); Chenu  _et al._ (2017); Cleve and Wang (2016) for Hamiltonian or
Lindbladian simulation problems. For gaussian bosonic environments Budini
(2000); Megier  _et al._ (2018), there are two prominent approaches to
identifying these environment modes. First is to use the pseudomode theory
Tamascelli  _et al._ (2018); Pleasance  _et al._ (2020); Mazzola  _et al._
(2009); Dalton  _et al._ (2001); Garraway and Dalton (2006); Dalton  _et al._
(2012), wherein the non-Markovian environment is approximated by a finite
number of discrete bosonic modes individually coupling to Markovian reservoirs
yielding an effective Lindbladian simulation problem. The second method is to
use Wilson’s start-to-delta transformation Chin  _et al._ (2010); Woods  _et
al._ (2014), which uses the Lanczos iteration to identify discrete bosonic
modes that approximate the environment and map the problem of computing non-
Markovian quantum dynamics to a Hamiltonian simulation problem.

While previous works have outlined a systematic construction of these
approximations, a rigorous theoretical study of their convergence properties
is lacking. In this paper, we develop a convergence theory of both the
pseudomode approximation as well as Wilson’s star-to-delta transformation when
applied to non-Markovian open quantum systems. We show that, for a large class
of non-Markovian systems with a Gaussian bosonic environment, the finite-time
dynamics of the quantum system can always be described, to a specified
accuracy, by either of these methods with a sufficiently large number of
discrete bosonic modes. Furthermore, we argue that for most non-Markovian
systems encountered in practice, the approximation error incurred falls off
polynomially with the number of environment modes. Our theoretical results
lend rigor to using these methods to approximate non-Markovian quantum
dynamics on classical, or quantum computers.

In order to perform this convergence study, there are several theoretical
challenges that our work resolves. First, since it is intuitively expected
that any approximation to a non-Markovian environment with only a finite
number of modes can only describe it within a finite bandwidth, we rigorously
argue that high energy components of the environment’s quantum state can be
neglected while studying finite-time dynamics of the quantum system. While
this would be relatively trivial to show for problems where the system-
environment couplings vanish at high energies, a number of models considered
in practice have non-vanishing couplings with arbitrarily high energy modes in
the environment Calajó  _et al._ (2019); Fan  _et al._ (2010) and it is the
magnitude of these couplings that limit the energies of the environment modes
that are excited. Furthermore, there is also a possibility of an ultraviolet
divergence, such as a diverging Lamb shift Barut and Kraus (1983) in the
system’s resonances, for several environment models which makes a rigorous
convergence study of the pseudomode approximation difficult. We resolve these
issues by identifying a set of physically motivated sufficient mathematical
conditions on the system-environment dynamics (assumptions 1 and 2 in section
II.2) that allow for neglecting the high energy components of the
environment’s quantum state (theorem 1 in section II.3). Building on this
result, we provide a bound on the error incurred in approximating the dynamics
of a non-Markovian quantum system, with system-environment couplings that
donot necessarily vanish at high energies, through the dynamics of a non-
Markovian quantum system wherein the system-environment couplings vanish at
high energies.

Next, we rigorously study the pseudomode approximation, and provide explicit
upper bounds on the error incurred in approximating a non-Markovian model with
pseudomodes. We show that this approximation becomes increasingly accurate
with the number of pseudomodes, thereby providing a convergence guarantee for
the pseudomode approximation (theorem 2 in section II.3). Furthermore, we
argue that for most models occurring in practice the rate of convergence is
polynomial with the number of pseudomodes used i.e. the approximation error
falls off polynomially with the number of pseudomodes. We derive a similar set
of results providing a convergence guarantee as well as an upper bound on the
approximation error for Wilson’s star-to-delta transformation (theorem 3 in
section II.3).

##  II Setup and Summary of results

###  II.1 Notation and mathematical prelimnaries

_Function spaces and analysis_ : For $v\in\mathbb{R}^{N}$ and
$u\in\mathbb{R}^{M}$, we will denote by $\\{v,u\\}$ a vector in
$\mathbb{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 $\mathbb{R}^{N}$. For a measure space
$(\Omega,\mathfrak{F},\mu)$, given a measurable map $f:\Omega\to\mathbb{R}$,
we define its essential supremum via
$\text{esssup}_{x\in\Omega}f(x)=\inf\\{a\in\mathbb{R}\ |\
\mu(f^{-1}([a,\infty]))=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
$\textnormal{P}:\Omega\to\\{\text{True},\text{False}\\}$ is said to be True
_almost everywhere_ in $\Omega$ if the $\mu\big{(}\\{x\ |\
\textnormal{P}(x)\text{ is False}\\}\big{)}=0$. This will be denoted as
$\forall x\in_{a.e.}\Omega:\ \text{P}(x)\text{ is True}$.

We will denote by $\text{L}^{2}(\mathbb{R}^{D})$ the set of complex-valued
square integrable functions over $\mathbb{R}^{D}$. We will denote by
$\text{B}(\mathbb{R}^{D})$ the set of bounded functions over $\mathbb{R}^{D}$.
For $f\in\text{B}(\mathbb{R}^{D})$, we denote by $\left\|{f}\right\|_{\infty}$
the maximum magnitude of the function i.e.
$\left\|{f}\right\|_{\infty}=\sup_{x\in\mathbb{R}^{D}}|f(x)|$. For a closed
interval $I\subset\mathbb{R}$, we will denote by $\text{C}^{k}(I)$ the set of
functions $f:I\to\mathbb{C}$ which are $k-$continuously differentiable. In
particular, $\text{C}^{0}(I)$ is the set of continuous complex-valued
functions with domain $I$ and $\text{C}^{\infty}(I)$ is the space of smooth
complex-valued functions with domain $I$. We will also denote by
$\textrm{C}^{k}_{b}(I)$ the space of bounded $k-$times continuously
differentiable functions — on compact intervals $I$, this will coincide with
$\textrm{C}^{k}(I)$. Finally, for any $\alpha>0$,
$\text{rect}_{\alpha}:\mathbb{R}\to\mathbb{R}$ to be denote the rectangular
windowing function defined by

$\text{rect}_{\alpha}(x)=\begin{cases}1&\text{ if }|x|\leq\alpha\\\ 0&\text{
if }|x|>\alpha.\end{cases}$

For $a<b$, we will denote by $\text{AC}([a,b])$ the space of all functions
$f:[a,b]\to\mathbb{C}$ that are absolutely continuous. We recall that
$f\in\text{C}^{0}([a,b])$ is absolutely continuous if $\exists
g:[a,b]\to\mathbb{C}$ such that $g$ is integrable and

$\forall\ x\in[a,b]:\ f(x)=f(a)+\int_{a}^{x}g(x^{\prime})dx^{\prime}.$

An immediate consequence of this definition is that the function $f$ has a
derivative, which we denote by $\partial_{x}f$, almost everywhere on $[a,b]$.
Furthermore, we can define a norm
$\left\|{\cdot}\right\|^{\mathcal{S}}_{[a,b]}:\text{AC}([a,b])\to\mathbb{R}^{+}$
such that $\forall f\in\text{AC}([a,b])$

$\displaystyle\left\|{f}\right\|^{\mathcal{S}}_{[a,b]}=\sup_{x\in[a,b]}|f(x)|+(b-a)\underset{x\in[a,b]}{\text{
esssup }}\big{|}\partial_{x}f(x)\big{|}.$ (1)

A map $f:[a,b]^{2}\to\mathbb{C}$ will be called symmetric if $\forall
s_{1},s_{2}\in[a,b]:\ f(s_{1},s_{2})=f^{*}(s_{2},s_{1})$. We will denote by
$\text{AC}_{\text{sym}}([a,b]^{2})$ the space of symmetric, continuous
functions $f:[a,b]^{2}\to\mathbb{C}$ which satisfy the following two
conditions

  1. 1.

$f$ is absolutely continuous in either of its arguments i.e. $\forall
y\in[a,b],$ the map $\varphi_{y}:[a,b]\to\mathbb{C}$ defined by
$\varphi_{y}(x)=f(x,y)$ is absolutely continuous.

  2. 2.

For $y\in_{a.e.}[a,b]$, $\partial_{y}f(x,y)$ exists $\forall x\in[a,b]$ and
the map $\varphi_{y}^{d}:[a,b]\to\mathbb{C}$ defined by
$\varphi_{y}^{d}(x)=\partial_{y}f(x,y)$ is absolutely continuous.

We note that the partial derivatives of $f$ with respect to either and both of
its arguments, $\partial_{x}f$, $\partial_{y}f$, $\partial^{2}_{x,y}f$, exists
almost everywhere on $[a,b]^{2}$. Finally, we define a norm
$\left\|{\cdot}\right\|_{[a,b]^{2}}^{\mathcal{S}}:\text{AC}_{\text{sym}}([a,b]^{2})\to\mathbb{R}^{+}$
such that $\forall f\in\text{AC}_{\text{sym}}([a,b]^{2})$

$\displaystyle\left\|{f}\right\|^{\mathcal{S}}_{[a,b]^{2}}$
$\displaystyle=\sup_{x,y\in[a,b]}|f(x,y)|+$
$\displaystyle(b-a)\bigg{(}\underset{x,y\in[a,b]}{\text{esssup}}\big{|}\partial_{x}f(x,y)\big{|}+\underset{{x,y\in[a,b]}}{\text{esssup}}\big{|}\partial_{y}f(x,y)\big{|}\bigg{)}+$
$\displaystyle(b-a)^{2}\underset{{x,y\in[a,b]}}{\text{esssup}}\big{|}\partial^{2}_{x,y}f(x,y)\big{|}.$
(2)

A symmetric map $K$ is also positive if

$\forall f\in\text{C}^{0}([a,b]):\
\int_{s_{1},s_{2}=a}^{b}K(s_{1},s_{2})f^{*}(s_{1})f(s_{2})ds_{1}ds_{2}\geq 0.$

_Norms_ : For a Hilbert space $\mathcal{H}$, the Banach space of bounded
operators over $\mathcal{H}$ will be denoted by $\mathfrak{L}(\mathcal{H})$.
For elements of a Hilbert space $\mathcal{H}$, $\ket{\psi}\in\mathcal{H}$,
$\left\|{\ket{\psi}}\right\|$ will denote the norm induced by
$\mathcal{H}^{\prime}$s inner product. For a bounded operator
$O:\mathcal{H}\to\mathcal{H}$, we will denote by $\left\|{O}\right\|$ the
operator norm induced by the norm on $\mathcal{H}$ i.e.
$\left\|{O}\right\|=\sup_{\ket{\psi}\in\mathcal{H}}\left\|{O\ket{\psi}}\right\|/\left\|{\ket{\psi}}\right\|$.
For a bounded and trace-class operator $O:\mathcal{H}\to\mathcal{H}$, its
trace-norm is given by
$\left\|{O}\right\|_{\text{tr}}=\text{tr}[\sqrt{O^{\dagger}O}]$. The trace-
norm will be used to quantify distance between two density operators. It is
useful to note that given two states
$\ket{\psi_{1}},\ket{\psi_{2}}\in\mathcal{H}$ with unit norm, then the
operator $\ket{\psi_{1}}\bra{\psi_{1}}-\ket{\psi_{2}}\bra{\psi_{2}}$ is
bounded and trace-class and its trace-norm satisfies

$\displaystyle\left\|{\ket{\psi_{1}}\bra{\psi_{1}}-\ket{\psi_{2}}\bra{\psi_{2}}}\right\|_{\text{tr}}\leq\sqrt{2}\left\|{\ket{\psi_{1}}-\ket{\psi_{2}}}\right\|.$
(3)

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:[0,\infty)\to[0,\infty)$ of
parameter $h$ then the notation $f=\mathcal{O}(g)\text{ as }h\to
h_{0}\implies\exists c>0:f(h)\leq cg(h)$ as $h\to h_{0}$. The notation
$f=\Theta(g)$ as $h\to h_{0}\implies\exists c_{-},c_{+}>0:c_{-}g(h)\leq
f(h)\leq c_{+}g(h)$ as $h\to h_{0}$. Furthermore,
$f=\mathcal{O}(g_{1})+\mathcal{O}(g_{2})$ as $h\to h_{0}$ if $\exists
f_{1},f_{2}$ such that $f=f_{1}+f_{2}$ and $f_{1}=\mathcal{O}(g_{1})$ as $h\to
h_{0}$ and $f_{2}=\mathcal{O}(g_{2})$ as $h\to h_{0}$.

###  II.2 Model and assumptions

We consider an open quantum system model, where a finite-dimensional quantum
system with hilbert space $\mathcal{H}_{S}=\mathbb{C}^{D}$ (referred to as a
‘local system’) interacting with a gaussian environment whose Hilbert space,
$\mathcal{H}_{B}$ is assumed to be a fock space over
$\text{L}^{2}(\mathbb{R})$. We denote by $a_{\omega}$ the annihilation
operator for this Fock space and consider hamiltonians of the form

$\displaystyle H(t)=$ $\displaystyle H_{S}(t)+\int_{-\infty}^{\infty}\omega
a^{\dagger}_{\omega}a_{\omega}d\omega+$
$\displaystyle\int_{-\infty}^{\infty}\bigg{(}v(\omega)a_{\omega}L^{\dagger}+v^{*}(\omega)a^{\dagger}_{\omega}L\bigg{)}d\omega,$
(4)

where $L:\mathcal{H}_{S}\to\mathcal{H}_{S}$ is the operator describing the
coupling of the system with the environment, and
$v\in\textrm{C}_{b}^{\infty}(\mathbb{R})$ is the frequency dependent coupling
function between the environment and the system. We point out that the
Hamiltonian in Eq. II.2 is only a provably (essentially) self-adjoint operator
if $v\in\text{L}^{2}(\mathbb{R})$ — however, several problems in quantum
optics (e.g. systems with point coupling) are described by $v$ that are
tempered distributions. To the best of our knowledge, there exists no
mathematically rigorous delineation of coupling functions $v$ that are not
square integrable but still result in a self-adjoint Hamiltonian, although
several such models are used in practice. A rigorous description of such
systems is likely to be constructed as some limit of systems described by
Hamiltonians of the form of Eq. II.2 with $v$ being in
$\text{L}^{2}(\mathbb{R})$, but is outside the scope of this paper and left
for future work. In this paper, we will assume without proof that the
Hamiltonian being considered has well defined dynamics.

Throughout this paper, we will be interested in approximating the dynamics of
the reduced state of the local system — we restrict ourselves to the case
where the environment is initially in a vacuum state, although an extension to
an initially excited environment state is easily possible. More explicitly,
starting from an initial state
$\ket{\psi(0)}=\ket{\sigma}\otimes\ket{\text{vac}}$, where
$\ket{\sigma}\in\mathcal{H}_{S}$ and $\ket{\text{vac}}\in\mathcal{H}_{B}$ is
the vacuum state of the environment, we are interested in computing
$\rho(t)=\text{Tr}_{B}[\ket{\psi(t)}\bra{\psi(t)}]$ where
$\ket{\psi(t)}=\mathcal{T}\exp(-i\int_{0}^{t}H(s)ds)\ket{\psi(0)}$.

An immediate issue with performing a rigorous theoretical analysis on the
convergence of pseudo-mode theory in approximating $\rho(t)$ is that
pseudomodes can only approximate the environment within a finite frequency
window. For most physically relevant models, we expect that the local system
cannot excite arbitrarily high frequencies in the environment in finite amount
of time even if it has a non-vanishing coupling with all frequencies in the
environment. Below we outline two physically motivated assumptions that we
make on the model under consideration, which can also be proved for a number
of physically interesting problems, that are sufficient conditions to allow us
to neglect high frequencies in the environment. While the first assumption is
on the coupling function, and can consequently be easily verified, the second
assumption is on the joint dynamics of the local system and the environment.

###### Assumption 1

The coupling function $v$ is smooth and bounded i.e.
$v\in\textnormal{C}_{b}^{\infty}(\mathbb{R})$. Furthermore, $\forall\ t\geq 0$
and $K\in\textnormal{AC}_{\textnormal{sym}}([0,t]^{2})$ such that $K$ is
positive,

$\displaystyle\lim_{\Omega\to\infty}\int_{\omega_{c}\leq|\omega_{1}|,|\omega_{2}|\leq\Omega}|v(\omega_{1})|^{2}|v(\omega_{2})|^{2}\bigg{(}\int_{s_{1},s_{2}=0}^{t}K(s_{1},s_{2})e^{i\omega_{2}(t-s_{2})}e^{-i\omega_{1}(t-s_{1})}ds\bigg{)}d\omega_{1}d\omega_{2}\leq
V(\omega_{c},t)\left\|{K}\right\|^{\mathcal{S}}_{[0,t]^{2}},$

where $V(\omega_{c},t)$ is a function that vanishes as $\omega_{c}\to\infty\
\forall t\geq 0$.

From a physical standpoint, this assumption ensures that the dynamics of the
model described by Eq. II.2 do not suffer from an ultraviolet divergence
arising due to the environment being able to support arbitrarily high
frequencies. We point out that it is unclear if this assumption, by itself, is
a necessary or sufficient condition to ensure that the dynamics generated by
the Hamiltonian in Eq. II.2 are well defined, and we leave this investigation
for future work. However, a number of models commonly considered in practice
do satisfy this condition:

  1. 1.

_Environments with square-integrable coupling functions:_ environments for
which the coupling function $v$, in addition to being smooth and bounded, is
square integrable i.e.
$v\in\textrm{C}_{b}^{\infty}(\mathbb{R})\cap\textrm{L}^{2}(\mathbb{R})$
satisfy assumption 1 (proposition 1 in appendix A). 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 Garraway (1997).

  2. 2.

_Markovian environments:_ Environments described by frequency-independent
coupling functions, i.e. $v(\omega)=v_{0}$, which are known to effectively
generate Markovian dynamics in the local system Fan  _et al._ (2010), also
satisfy assumption 1 (proposition 2 in appendix A).

  3. 3.

_Environments modelling retardation effects:_ environments described by
frequency-dependent coupling constants of the form
$v(\omega)=\sum_{i=1}^{M}v_{i}e^{i\omega\tau_{i}}$ are commonly used to model
retardation effects Cilluffo  _et al._ (2020). These environments also satisfy
the conditions of assumption 1 (proposition 2 in appendix A).

The second assumption that we introduce is related to the joint dynamics of
the local system and the environment. To state this assumption mathematically,
we introduce the $N-$point Green’s function of the localized system evaluated
at $s\in[0,t]^{N}$, $G_{N}(s;t)\in\mathfrak{L}(\mathcal{H}_{S})$ via

$\displaystyle
G_{N}(s;t)=\bra{\text{vac}}U(t,0)\mathcal{T}\bigg{[}\prod_{i=1}^{N}L(s_{i})\bigg{]}\ket{\text{vac}},$
(5)

where $U(\tau_{1},\tau_{2})$ is the propagator corresponding to the
Hamiltonian in Eq. II.2 and $L(\tau)=U(0,\tau)LU(\tau,0)$ is the operator $L$
in the Heisenberg picture. It can be immediately noted that the Green’s
function is a bounded operator, and satisfies the estimate

$\displaystyle\forall
s\in[0,t]^{N}:\left\|{G_{N}(s;t)}\right\|_{2}\leq\left\|{L}\right\|_{2}^{N}.$
(6)

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. More explicitly, consider a setting where the system is initialized
in a state $\ket{\sigma}$ and the environment is in its vacuum state
$\ket{\text{vac}}$, then the state of the system at time $t$ can be expressed
(see appendix B for a derivation):

$\displaystyle\ket{\psi(t)}=\sum_{N=0}^{\infty}\ket{\psi^{N}(t)},$ (7a) where
$\displaystyle\ket{\psi^{N}(t)}=\frac{1}{N!}\int_{\omega\in\mathbb{R}}F_{N}(\omega;t)\bigg{[}\prod_{i=1}^{N}v^{*}(\omega_{i})a_{\omega_{i}}^{\dagger}\bigg{]}\ket{\sigma}\otimes\ket{\text{vac}}d\omega,$
(7b) and $\forall\omega\in\mathbb{R}^{N}$ $\displaystyle
F_{N}(\omega;t)=\int_{s\in[0,t]^{N}}G_{N}(s;t)e^{-i\omega\cdot(t-s)}ds.$ (7c)

Furthermore, the bound on the $N-$point Green’s function (Eq. 6) implies the
following bound on the norm of $\ket{\psi^{N}(t)}\ \forall t\geq 0,N\geq 0$:

$\displaystyle\left\|{\ket{\psi^{N}(t)}}\right\|^{2}\leq\frac{\left\|{v}\right\|_{\infty}^{2N}\left\|{L}\right\|_{2}^{2N}(2\pi
t)^{N}}{N!}.$ (8)

###### Assumption 2

$\forall t\geq 0,s\in[0,t]^{N-1}$, the map
$G_{N}(\\{\cdot,s\\};t):[0,t]\to\mathfrak{L}(\mathcal{H}_{S})$ is absolutely
continuous and $\exists\gamma(t)>0$ such that

$\displaystyle\underset{s_{0}\in[0,t]}{\textnormal{esssup}}\left\|{\partial_{s_{0}}G_{N}(\\{s_{0},s\\};t)}\right\|_{2}\leq\left\|{L}\right\|_{2}^{N}\gamma(t).$

This assumption can be interpreted as a bound on the rate at which the local
system can emit an excitation into the environment, despite its non-negligible
couplings with arbitrarily high frequencies in the environment modes. While we
expect this assumption to be true for all physically reasonable open quantum
system models, this assumption can be proved for two cases:

  1. 1.

_A Markovian environment_ , i.e. environments with a frequency independent
coupling constant: $v(\omega)=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 appendix C).

  2. 2.

_A environment with a square integrable coupling constant:_ For environments
with a bounded and square integrable coupling constant $v$, assumption 2 can
again be rigorously proven (proposition 4 in appendix C).

###  II.3 Summary of results

With these two assumptions, we can now prove the convergence of the pseudomode
theory Tamascelli  _et al._ (2018) and Wilson’s start-to-delta transformation
Chin  _et al._ (2010) for simulating non-Markovian quantum systems. Our first
result rigorously shows that in a finite amount of time, the localized system
cannot excite arbitrarily high frequencies in the environment. More precisely,

###### Theorem 1

Suppose $v\in\textnormal{C}_{b}^{\infty}(\mathbb{R})$ is a coupling function
such that assumptions 1 and 2 are satisfied. Denoting by $\rho(t)$ the reduced
density matrix of the local system at time $t$ when an initial state
$\ket{\sigma}\otimes\ket{\text{vac}}$ is evolved under the Hamiltonian in Eq.
II.2 and by ${\rho}_{\omega_{c}}(t)$ the reduced density matrix of the local
system at time $t$ when the same initial state is evolved under the
Hamiltonian

$\displaystyle H_{\omega_{c}}(t)=$ $\displaystyle
H_{S}(t)+\int_{-\infty}^{\infty}\omega a^{\dagger}_{\omega}a_{\omega}d\omega+$
$\displaystyle\int_{-\infty}^{\infty}\bigg{(}v(\omega)a_{\omega}L^{\dagger}+v^{*}(\omega)a^{\dagger}_{\omega}L\bigg{)}\textnormal{rect}_{\omega_{c}}(\omega)d\omega,$
(9)

then

$\displaystyle\left\|{\rho(t)-{\rho}_{\omega_{c}}(t)}\right\|_{\textnormal{tr}}\leq\varepsilon(\omega_{c},t),$

where $\varepsilon(\omega_{c},t)$ is the cutoff error given by

$\displaystyle\varepsilon(\omega_{c},t)=\frac{f_{1}(t)}{\sqrt{\omega_{c}}}+\int_{0}^{t}f_{2}(\tau)\sqrt{V(\omega_{c},\tau)}d\tau.$
(10)

Here $V(\omega_{c},t)$ is defined in assumption 1 and

$\displaystyle
f_{1}(\tau)=\sqrt{2}\left\|{v}\right\|_{\infty}\left\|{L}\right\|_{2}(2+\gamma(\tau)\tau)e^{\left\|{v}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}\tau/2},$
$\displaystyle
f_{2}(\tau)=\sqrt{2}\left\|{L}\right\|_{2}^{2}\tau(1+\gamma(\tau)\tau)e^{\left\|{v}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}\tau/2},$

where $\gamma(t)$ is introduced in assumption 2.

Consequently, the environment can be approximated by its restriction within a
finite frequency window. While we will make use of this result in proving the
convergence of the pseudomode approximation below, we point out that this
result also lends theoretical rigor to a number of other numerical methods
used in practice which either explicitly introduce a frequency cutoff on the
environment modes Chin  _et al._ (2010); Woods  _et al._ (2014), or implicitly
introduce a frequency cutoff by discretizing the environment modes in
position-domain Fang (2019).

We point out that the bounds that we provide in theorem 1 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 excitations into the
environment, and the only bound on the probability of emission of $N$
particles that can be derived generally, as given by Eq. 8, is poissonian in
$N$ with a parameter that grows with $t$. In practice, we expect the errors to
not grow exponentially with $t$ and the scaling of this error with $t$ to
depend on the number of particles that the local system is emitting into the
environment.

Next, we consider a concrete approximation method based on pseudo-mode theory.
Here, for a given coupling function $v$, we approximate $|v|^{2}$ by a sum of
lorentzians within a finite frequency window, and each Lorentzian can be
described as an independent discrete bosonic mode coupling to a Markovian
environment. It thus approximates the dynamics of the local system when
evolved according to the Hamiltonian in Eq. II.2 by the dynamics of a larger
Markovian system. More precisely, for a pseudomode description with $M$
bosonic modes, we consider an effective system whose Hilbert space is
$\mathcal{H}_{S}\otimes\mathcal{H}_{\text{aux}}$ where
$\mathcal{H}_{\text{aux}}$ is a Fock space over $\mathbb{C}^{M}$, with a
Hamiltonian $\hat{H}(t)$ given by

$\displaystyle\hat{H}(t)=H_{S}(t)+\sum_{i=1}^{M}\omega_{i}a_{i}^{\dagger}a_{i}+\sum_{i=1}^{M}g_{i}(a_{i}L^{\dagger}+a_{i}^{\dagger}L).$

Here $\omega_{i}$ and $g_{i}$ for $i\in\\{1,2\dots M\\}$ are parameters of the
pseudomode description. The time-evolution of the effective system is governed
by a Lindbladian master equation for its density matrix $R(t)$:

$\displaystyle\frac{dR(t)}{dt}=i[\hat{H}(t),R(t)]+\sum_{i=1}^{M}\frac{\kappa_{i}}{2}\big{(}2a_{i}R(t)a_{i}^{\dagger}-\\{a_{i}^{\dagger}a_{i},R(t)\\}\big{)}.$

where $a_{i}$ are the annihilation operators corresponding to the introduced
bosonic modes and $\kappa_{i}$ for $i\in\\{1,2\dots M\\}$ are their decay
rates which are also parameters of the pseudomode description. This pseudomode
description effectively computes the dynamics of the local system coupling to
a bath with a coupling constant $\hat{v}$ where Tamascelli  _et al._ (2018)

$|\hat{v}(\omega)|^{2}=\sum_{i=1}^{M}\frac{\kappa_{i}}{2\pi}\frac{g_{i}^{2}}{(\omega-\omega_{i})^{2}+\kappa_{i}^{2}/4}.$

To setup the pseudomode description for a coupling constant $v$, we choose the
parameters $\omega_{i},g_{i},\kappa_{i}$ for $i\in\\{1,2\dots M\\}$ to best
approximate $|\hat{v}|^{2}$ with $|v|^{2}$ within a sufficiently large
frequency window. The reduced state of the local system at time $t$, obtained
on evolving an initial state $\ket{\sigma}\otimes\ket{\text{vac}},\
\ket{\sigma}\in\mathcal{H}_{S}$, can then be computed by evolving an initial
density matrix $R(0)=\ket{\sigma}\bra{\sigma}\otimes(\ket{0}\bra{0})^{\otimes
M}$ using this master equation and finally tracing over the auxillary system.

We show that, for coupling functions that satisfy the two assumptions laid out
in the previous subsection, pseudomode approximation of the non-Markovian
system for computing its finite-time dynamics is guaranteed to converge for
sufficiently large number of pseudomodes. Furthermore, we also characterize
the rate of convergence i.e. how fast does the error incurred in the
pseudomode approximation decrease with the number of pseudomodes used in the
approximation. The convergence of the pseudomode approximation will, in
general, depend on the details of the coupling function $v$ — in particular,
on how fast its derivative $v^{\prime}(\omega)$ can grow with $\omega$, as
well as how the cutoff error introduced in theorem 1 decreases with the cutoff
frequency. For typical coupling functions encountered in practice, we expect a
polynomial growth in the derivative of the coupling function with frequency
and a polynomial fall off in the cutoff error with the cutoff frequency. Under
these assumptions, as detailed in theorem below, we can show that the error
incurred in the pseudomode approximation falls over polynomially in the number
of pseudomodes.

######  Theorem 2 (Pseudomode approximation)

Suppose $v\in\textnormal{C}^{\infty}_{b}(\mathbb{R})$ is a coupling function
such that assumptions 1 and 2 are satisfied and let $\rho(t)$ be the reduced
state of the local system after evolving an initial state
$\ket{\sigma}\otimes\ket{\textnormal{vac}}$ using the Hamiltonian in Eq. II.2.
Then, there exists a pseudomode description of the non-Markovian system with
$M$ bosonic modes which provides an approximation $\hat{\rho}(t)$ to the
reduced state of the local system such that
$\left\|{\rho(t)-\hat{\rho}(t)}\right\|_{\textnormal{tr}}\to 0$ as
$M\to\infty$. Furthermore, if
$|v^{\prime}(\omega)|=O(\textnormal{poly}(\omega))$ and the cutoff error
$\varepsilon(\omega_{c},t)=\mathcal{O}(\textnormal{poly}(t)\exp(\mathcal{O}(t))\textnormal{poly}(\omega_{c}^{-1}))$,
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
$\mathcal{O}(\textnormal{poly}(t)\exp(\mathcal{O}(t))\textnormal{poly}(M^{-1}))$.

Finally, we consider the Wilson star-to-delta transformation, which utilizes a
Lanczos recurrence Van Dooren (1995) to approximate the non-Markovian
environment. The approximation first introduces a frequency cutoff
$\omega_{c}$ on the environment modes. Then, starting from an initial mode
$A_{0}\propto\int_{-\infty}^{\omega_{c}}v(\omega)a_{\omega}d\omega$, the
Lanczos iteration with respect to the environment Hamiltonian is used to
generate modes $A_{1},A_{2},A_{3}\dots$ with which the environment can be
modelled by a 1D nearest neighbour model with non-uniform coupling strengths

$H^{E}_{\omega_{c}}=\sum_{i=0}^{\infty}\omega_{i}A_{i}^{\dagger}A_{i}+\sum_{i=0}^{\infty}g_{i}(A_{i}A_{i+1}^{\dagger}+A_{i}^{\dagger}A_{i+1}),$

where the parameters $\omega_{i}$ and $g_{i}$ for $i\in\\{0,1,2\dots\\}$ are
obtained from the Lanczos iteration (see section III.3 for details).
Furthermore, accounting for the frequency cutoff, the interaction Hamiltonian
in between the local system and the environment is
$\propto(A_{0}L^{\dagger}+LA_{0}^{\dagger})$. If all the (infinite number of)
modes $A_{0},A_{1},A_{2}\dots$ are accounted for, then within the frequency
cutoff, the dynamics of the non-Markovian system are captured exactly.
However, if only a finite number of modes are accounted for (as in a practical
numerical simulation), then an error will be incurred. We show that the, given
assumptions 1 and 2 hold, the Wilson start-to-delta transformation will
converge for sufficiently large number of modes and that the convergence rate
is polynomial in the number of modes.

######  Theorem 3 (Wilson star-to-delta transformation)

Suppose $v\in\textnormal{C}^{\infty}_{b}(\mathbb{R})$ is a coupling function
such that assumptions 1 and 2 are satisfied and let $\rho(t)$ be the reduced
state of the local system after evolving an initial state
$\ket{\sigma}\otimes\ket{\textnormal{vac}}$ using the Hamiltonian in Eq. II.2.
Then, using Wilson’s star-to-delta transformation with $M$ modes provides an
approximation $\hat{\rho}(t)$ to the reduced state of the local system such
that $\left\|{\rho(t)-\hat{\rho}(t)}\right\|_{\textnormal{tr}}\to 0$ as
$M\to\infty$. Furthermore, if the cutoff error
$\varepsilon(\omega_{c},t)=\mathcal{O}(\textnormal{poly}(t)\exp(\mathcal{O}(t))\textnormal{poly}(\omega_{c}^{-1}))$,
then the trace-norm error in approximating the reduced local system state at
time $t$ scales as
$\mathcal{O}(\textnormal{poly}(t)\exp(\mathcal{O}(t))\textnormal{poly}(M^{-1}))$.

Our proof of this convergence result relies on the connection between the
Wilson’s star-to-delta transformation and orthogonal polynomials that has been
previously established Chin  _et al._ (2010); de Vega  _et al._ (2015).

The remainder of this paper contains detailed proofs of theorems 1–3. We only
prove the most important results in the main text, and relegate longer
calculations to the appendices.

##  III Detailed Proofs

###  III.1 Introducing a frequency cut-off on the environment

We first show rigorously that, subject to the assumptions 1 and 2 stated in
Section II.2, a frequency cutoff can be introduced in the environment without
significantly impacting the dynamics of the local system. For a given
frequency cutoff $\omega_{c}>0$, we introduce an operator,
$\mathcal{P}_{\omega_{c}}:\mathcal{H}_{S}\otimes\mathcal{H}_{B}\to\mathcal{H}_{S}\otimes\mathcal{H}_{B}$,
that projects the joint system and environment state onto a state with no
environment excitations at frequencies $|\omega|\geq\omega_{c}$ — more
rigorously, $\forall N\geq
0,\sigma\in\mathcal{H}_{S},\psi\in\textrm{L}^{2}(\mathbb{R}^{N})$

$\displaystyle\mathcal{P}_{\omega_{c}}\int_{\omega\in\mathbb{R}^{N}}\psi(\omega)\bigg{[}\prod_{i=1}^{N}a_{\omega_{i}}^{\dagger}\bigg{]}\ket{\sigma}\otimes\ket{\text{vac}}d\omega$
$\displaystyle\qquad=\int_{\omega\in[-\omega_{c},\omega_{c}]^{N}}\psi(\omega)\bigg{[}\prod_{i=1}^{N}a_{\omega_{i}}^{\dagger}\bigg{]}\ket{\sigma}\otimes\ket{\text{vac}}d\omega.$
(11)

Furthermore, we define
$\mathcal{Q}_{\omega_{c}}=\text{I}-\mathcal{P}_{\omega_{c}}$. It is
straightforward to see that the Hamiltonians defined in Eqs. II.2 and 1
satisfy
$H_{\omega_{c}}(t)=\mathcal{P}_{\omega_{c}}H(t)\mathcal{P}_{\omega_{c}}$. We
now study the error introduced by evolving an initial state under the
Hamiltonian $H_{\omega_{c}}(t)$ instead of the Hamiltonian $H(t)$. This is
made explicit in the following lemma.

###### Lemma 1

Given an initial state $\ket{\psi_{0}}$ that satisfies
$\mathcal{P}_{\omega_{c}}\ket{\psi_{0}}=\ket{\psi_{0}}$, and denoting by
$U(t_{1},t_{2})$ and $U_{\omega_{c}}(t_{1},t_{2})$ the propagators
corresponding to $H(t)$ and $H_{\omega_{c}}(t)$ respectively, the states
$\ket{\psi(t)}=U(t,0)\ket{\psi_{0}}$ and
$\ket{\psi_{\omega_{c}}(t)}=U_{\omega_{c}}(t,0)\ket{\psi_{0}}$ satisfy

$\displaystyle\left\|{\ket{\psi(t)}-\ket{\psi_{\omega_{c}}(t)}}\right\|\leq$
$\displaystyle\qquad\left\|{\mathcal{Q}_{\omega_{c}}\ket{\psi(t)}}\right\|+\int_{0}^{t}\left\|{\mathcal{P}_{\omega_{c}}H(s)\mathcal{Q}_{\omega_{c}}\ket{\psi(s)}}\right\|ds$

Proof: Note that
$\ket{\psi(t)}=\mathcal{P}_{\omega_{c}}\ket{\psi(t)}+\mathcal{Q}_{\omega_{c}}\ket{\psi(t)}$.
From the Schroedinger’s equation for $\ket{\psi(t)}$, it then follows that

$\displaystyle i\frac{d}{dt}\mathcal{P}_{\omega_{c}}\ket{\psi(t)}$
$\displaystyle=\mathcal{P}_{\omega_{c}}H(t)\ket{\psi(t)},$
$\displaystyle=H_{\omega_{c}}(t)\mathcal{P}_{\omega_{c}}\ket{\psi(t)}+\mathcal{P}_{\omega_{c}}H(t)\mathcal{Q}_{\omega_{c}}\ket{\psi(t)},$

which can be integrated to obtain

$\displaystyle\mathcal{P}_{\omega_{c}}\ket{\psi(t)}=\ket{\psi_{\omega_{c}}(t)}-i\int_{0}^{t}U_{\omega_{c}}(t,s)\mathcal{P}_{\omega_{c}}H(s)\mathcal{Q}_{\omega_{c}}\ket{\psi(s)}ds.$

Consequently

$\displaystyle\ket{\psi(t)}-\ket{\psi_{\omega_{c}}(t)}=$
$\displaystyle\qquad\mathcal{Q}_{\omega_{c}}\ket{\psi(t)}-i\int_{0}^{t}U_{\omega_{c}}(t,s)\mathcal{P}_{\omega_{c}}H(s)\mathcal{Q}_{\omega_{c}}\ket{\psi(s)}ds,$

from which the statement of the lemma follows. $\square$

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 $\omega_{c}$ is made increasingly larger.

###### Lemma 2

If assumptions 1 and 2 are satisfied, it follows that $\forall t>0$,

$\displaystyle\left\|{\mathcal{Q}_{\omega_{c}}\ket{\psi(t)}}\right\|_{2}\leq\frac{\left\|{v}\right\|_{\infty}\left\|{L}\right\|_{2}}{\sqrt{\omega_{c}}}(2+\gamma(t)t)e^{\pi\left\|{v}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}t},$

where $\gamma(t)$ is introduced in assumption 2.

Proof: It follows from the definition of $\mathcal{Q}_{\omega_{c}}$ that

$\left\|{Q_{\omega_{c}}\ket{\psi(t)}}\right\|_{2}=\sqrt{\bra{\psi(t)}(\text{I}-\mathcal{P}_{\omega_{c}})\ket{\psi(t)}}.$

Using Eq. 7, we obtain

$\displaystyle\bra{\psi(t)}\psi(t)\rangle=$
$\displaystyle\sum_{N=0}^{\infty}\frac{1}{N!}\int_{\mathbb{R}^{N}}\bra{\sigma}F_{N}^{\dagger}(\omega;t)F_{N}(\omega;t)\ket{\sigma}\prod_{i=1}^{N}|v(\omega_{i})|^{2}d\omega$

and

$\displaystyle\bra{\psi(t)}\mathcal{P}_{\omega_{c}}\ket{\psi(t)}=$
$\displaystyle\sum_{N=0}^{\infty}\frac{1}{N!}\int_{[-\omega_{c},\omega_{c}]^{N}}\bra{\sigma}F_{N}^{\dagger}(\omega;t)F_{N}(\omega;t)\ket{\sigma}\prod_{i=1}^{N}|v(\omega_{i})|^{2}d\omega.$
(12)

Therefore,

$\displaystyle\bra{\psi(t)}\mathcal{Q}_{\omega_{c}}\ket{\psi(t)}\leq$
$\displaystyle\sum_{N=1}^{\infty}\frac{\left\|{v}\right\|_{\infty}^{2N}}{(N-1)!}\int_{\begin{subarray}{c}|\omega_{1}|\geq\omega_{c}\\\
\omega_{2}\dots\omega_{N}\in\mathbb{R}\end{subarray}}\bra{\sigma}F_{N}^{\dagger}(\omega;t)F_{N}(\omega;t)\ket{\sigma}d\omega=$
$\displaystyle\sum_{N=1}^{\infty}\frac{(2\pi\left\|{v}\right\|_{\infty}^{2})^{N}}{2\pi(N-1)!}\int_{\begin{subarray}{c}|\omega|\geq\omega_{c}\\\
s\in[0,t]^{N-1}\end{subarray}}\bra{\sigma}H_{N}^{\dagger}(\omega,s;t)H_{N}(\omega,s;t)\ket{\sigma}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
d\omega\ ds.$

where, for $N\geq 1$ and $s\in[0,t]^{N-1}$, $H_{N}(\omega,s;t)$ is given by

$\displaystyle
H_{N}(\omega,s;t)=\int_{s^{\prime}\in[0,t]}G_{N}([s^{\prime},s];t)e^{-i\omega(t-s^{\prime})}ds^{\prime}.$

Using integration by parts, together with assumption 2, we immediately obtain
that

$\displaystyle\left\|{H_{N}(\omega,s;t)}\right\|_{2}\leq\frac{(2+\gamma(t)t)}{\omega}\left\|{L}\right\|_{2}^{N},$

and therefore

$\displaystyle\bra{\psi(t)}\mathcal{Q}_{\omega_{c}}\ket{\psi(t)}\leq\frac{\left\|{v}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}}{\omega_{c}}(2+\gamma(t)t)^{2}e^{2\pi\left\|{v}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}t}.$

from which the lemma immediately follows $\square$.

###### Lemma 3

If assumptions 1 and 2 are satisfied, it follows that $\forall t>0$,

$\displaystyle\int_{0}^{t}\left\|{\mathcal{P}_{\omega_{c}}H(\tau)\mathcal{Q}_{\omega_{c}}\ket{\psi(\tau)}}\right\|_{2}d\tau\leq$
$\displaystyle\qquad\left\|{L}\right\|_{2}^{2}\int_{0}^{t}\tau(1+\gamma(\tau)\tau)\big{[}V(\omega_{c},\tau)\big{]}^{1/2}e^{\pi\left\|{v}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}\tau}d\tau,$

where $\gamma(t)$ is introduced in assumption 2.

Proof: Using the definitions of $\mathcal{P}_{\omega_{c}}$ and
$\mathcal{Q}_{\omega_{c}}$, it follows that

$\displaystyle\mathcal{P}_{\omega_{c}}H(\tau)\mathcal{Q}_{\omega_{c}}\ket{\psi(\tau)}=$
$\displaystyle\qquad
L^{\dagger}\sum_{N=1}^{\infty}\frac{1}{(N-1)!}\int_{\begin{subarray}{c}|\omega_{1}|\geq\omega_{c}\\\
\omega_{2}\dots\omega_{N}\in\mathbb{R}\end{subarray}}|v(\omega_{1})|^{2}F_{N}(\omega;\tau)\times$
$\displaystyle\qquad\qquad\qquad\qquad\bigg{[}\prod_{i=2}^{N}v^{*}(\omega_{i})a_{\omega_{i}}^{\dagger}\bigg{]}\ket{\sigma}\otimes\ket{\text{vac}}d\omega.$

Therefore,

$\displaystyle\left\|{\mathcal{P}_{\omega_{c}}H(\tau)\mathcal{Q}_{\omega_{c}}\ket{\psi(\tau)}}\right\|^{2}$
$\displaystyle\leq\sum_{N=1}^{\infty}\frac{(2\pi\left\|{v}\right\|_{\infty}^{2})^{N-1}}{(N-1)!}\int_{\begin{subarray}{c}s\in[0,t]^{N-1}\\\
|\omega|,|\omega^{\prime}|\geq\omega_{c}\end{subarray}}|v(\omega)|^{2}|v(\omega^{\prime})|^{2}\times$
$\displaystyle\qquad\qquad\bra{\sigma}H_{N}^{\dagger}(\omega^{\prime},s;\tau)LL^{\dagger}H_{N}(\omega,s;\tau)\ket{\sigma}d\omega\
d\omega^{\prime}ds.$

To proceed further, we make use of assumption 1. For $\tau>0,N\geq 1,\
s\in[0,t]^{N-1}$, we construct a positive kernel
$K\in\text{AC}_{\text{sym}}([0,\tau]^{2})$ from the $N-$point Green’s function
via
$K(t^{\prime},t^{\prime\prime})=\bra{\sigma}G_{N}^{\dagger}(\\{t^{\prime},s\\};\tau)LL^{\dagger}G_{N}(\\{t^{\prime\prime},s\\};\tau)\ket{\sigma}$.
An application of the assumption 1 using this kernel yields $\forall
s\in[0,t]^{N-1}$

$\displaystyle\int_{|\omega|,|\omega^{\prime}|\geq\omega_{c}}|v(\omega)|^{2}|v(\omega^{\prime})|^{2}\times$
$\displaystyle\qquad\qquad\bra{\sigma}H_{N}^{\dagger}(\omega^{\prime},s;\tau)LL^{\dagger}H_{N}(\omega,s;\tau)\ket{\sigma}d\omega^{\prime}d\omega\leq$
$\displaystyle\qquad
V(\omega_{c},\tau)\left\|{\bra{\sigma}G_{N}^{\dagger}(\\{\cdot,s\\};\tau)LL^{\dagger}G_{N}(\\{\cdot,s\\};\tau)\ket{\sigma}}\right\|_{[0,\tau]^{2}}^{\mathcal{S}},$

where $V(\omega_{c},t)$ is the function introduced in assumption 1. Using
assumption 2, we immediately obtain that $\forall s\in[0,\tau]^{N-1}$

$\displaystyle\left\|{\bra{\sigma}G_{N}^{\dagger}([\cdot,s];\tau)LL^{\dagger}G_{N}([\cdot,s];\tau)\ket{\sigma}}\right\|_{[0,t]^{2}}^{\mathcal{S}}\leq$
$\displaystyle\qquad\ \ \ \left\|{L}\right\|_{2}^{2N+2}(1+\gamma(t)t)^{2}$

It thus follows that

$\displaystyle\left\|{\mathcal{P}_{\omega_{c}}H(\tau)\mathcal{Q}_{\omega_{c}}\ket{\psi(\tau)}}\right\|\leq\left\|{L}\right\|_{2}^{2}(1+\gamma(\tau)\tau)e^{\pi\left\|{v}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}\tau}$

Using this bound, the lemma statement follows. $\square$

Proof of theorem 1: Let $\ket{\sigma}\in\mathcal{H}_{S}$ be an arbitrary local
system state, and denote by $\ket{\psi(t)}$ the joint state of the system and
the environment obtained on evolving $\ket{\sigma}\otimes\ket{\text{vac}}$
under the Hamiltonian $H(t)$. Similarly, we denote by
$\ket{\psi_{\omega_{c}}(t)}$ the joint state of the system and environment
obtained on evolving $\ket{\sigma}\otimes\ket{\text{vac}}$ under the
Hamiltonian
$H_{\omega_{c}}(t)=\mathcal{P}_{\omega_{c}}H(t)\mathcal{P}_{\omega_{c}}$. We
denote by $\rho(t)$ and $\rho_{\omega_{c}}(t)$ the reduced state of the local
system at time $t$ i.e.
$\rho(t)=\text{Tr}_{\mathcal{H}_{B}}[\ket{\psi(t)}\bra{\psi(t)}]$ and
$\rho_{\omega_{c}}(t)=\text{Tr}_{\mathcal{H}_{B}}[\ket{\psi(t)}\bra{\psi(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

$\displaystyle\left\|{\rho(t)-\rho_{\omega_{c}}(t)}\right\|_{\text{tr}}$
$\displaystyle\leq\left\|{\ket{\psi(t)}\bra{\psi(t)}-\ket{\psi_{\omega_{c}}(t)}\bra{\psi_{\omega_{c}}(t)}}\right\|_{\text{tr}}$
$\displaystyle\leq\sqrt{2}\left\|{\ket{\psi(t)}-\ket{\psi_{\omega_{c}}(t)}}\right\|.$

An application of lemmas 1, 2 and 3 allows us to bound
$\left\|{\ket{\psi(t)}-\ket{\psi_{\omega_{c}}(t)}}\right\|$ and obtain the
result in theorem $\square$.

###  III.2 Convergence of pseudomode approximation

Having established that, for the class of open quantum systems under
consideration, high frequencies of the environment modes can be safely
neglected, we next consider approximating the environment by pseudomodes
within a sufficiently large bandwidth. In order to rigorously show that this
is possible, 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 4

Let
$v_{1},v_{2}\in\textnormal{C}^{\infty}_{b}(\mathbb{R})\cap\textnormal{L}^{2}(\mathbb{R})$
be coupling functions and let $H_{1}(t),H_{2}(t)$ denote the Hamiltonians
corresponding to these coupling functions as given by Eq. II.2. Denoting by
$\ket{\psi_{i}(t)}$ as the joint state of the local system and environment
when an initial state $\ket{\sigma}\otimes\ket{\textnormal{vac}}$, with
$\ket{\sigma}\in\mathcal{H}_{S}$, is evolved under the Hamiltonian $H_{i}(t),\
i\in\\{1,2\\}$, then

$\displaystyle\left\|{\ket{\psi_{1}(t)}-\ket{\psi_{2}(t)}}\right\|\leq
g(t)d_{\omega_{c}}(v_{1},v_{2}),$ (13)

where

$\displaystyle
g(t)=\int_{0}^{t}\bigg{(}\sqrt{1+\left\|{v_{1}}\right\|_{\infty}^{2}\tau\left\|{L}\right\|_{2}^{2}e^{2\pi\left\|{v_{1}}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}\tau}}+$
$\displaystyle\left\|{v_{1}}\right\|_{\infty}\sqrt{\tau}\left\|{L}\right\|_{2}e^{\pi\left\|{v_{1}}\right\|^{2}_{\infty}\left\|{L}\right\|_{2}^{2}\tau}\bigg{)}d\tau$

and

$d(v_{1},v_{2})=\bigg{[}\int_{-\infty}^{\infty}|v_{1}(\omega)-v_{2}(\omega)|^{2}d\omega\bigg{]}^{1/2},$

quantifies a distance between the two coupling functions.

Proof: For $i\in\\{1,2\\}$, we denote by $U_{i}(\tau_{1},\tau_{2})$ the
propagator corresponding to the Hamiltonian $H_{i}$. Furthermore, let
$\ket{\psi_{i}(t)}=U_{i}(t,0)\ket{\sigma}\otimes\ket{\text{vac}}$, where
$\ket{\sigma}\in\mathcal{H}_{S}$. It can be seen that

$\left\|{\ket{\psi_{1}(t)}-\ket{\psi_{2}(t)}}\right\|=\left\|{\ket{\phi(t)}-\ket{\sigma}\otimes\ket{\text{vac}}}\right\|,$

where

$\ket{\phi(t)}=U_{2}^{\dagger}(t,0)\ket{\psi_{1}(t)}.$

An integral equation can easily be derived by differentiating this equation,
followed by integrating it from $0$ to $t$ to obtain

$\ket{\phi(t)}=\ket{\sigma}\otimes\ket{\text{vac}}+i\int_{0}^{t}U_{2}(t,\tau)(H_{2}(\tau)-H_{1}(\tau))\ket{\psi_{1}(\tau)}d\tau,$

from which we obtain the bound

$\left\|{\psi_{1}(t)-\psi_{2}(t)}\right\|\leq\int_{0}^{t}\left\|{(H_{2}(\tau)-H_{1}(\tau))\ket{\psi_{1}(\tau)}}\right\|d\tau.$

We note from the definition of the Hamiltonians $H_{i}(t)$ that
$H_{2}(t)-H_{1}(t)=\Delta^{+}L+L^{\dagger}\Delta^{-}$ where

$\displaystyle\Delta^{+}=\int_{-\infty}^{\infty}\big{(}v_{1}^{*}(\omega)-v_{2}^{*}(\omega)\big{)}a_{\omega}^{\dagger}d\omega,$
$\displaystyle\Delta^{-}=\int_{-\infty}^{\infty}\big{(}v_{1}(\omega)-v_{2}(\omega)\big{)}a_{\omega}d\omega.$

Therefore,

$\displaystyle\left\|{\psi_{1}(t)-\psi_{2}(t)}\right\|\leq$
$\displaystyle\qquad\left\|{L}\right\|_{2}\bigg{(}\int_{0}^{t}\left\|{\Delta^{+}\ket{\psi_{1}(\tau)}}\right\|d\tau+\int_{0}^{t}\left\|{\Delta^{-}\ket{\psi_{1}(\tau)}}\right\|d\tau\bigg{)}.$
(14)

Note that the operators $\Delta^{\pm}$ are unbounded, and consequently the
norms in the above equations cannot be estimated trivially by introducing the
operator norms of $\Delta^{\pm}$. However, $\Delta^{\pm}$ are bounded when
restricted to only the $N-$particle subspaces of the environment’s Hilbert
space. In particular, for a $N-$particle state $\ket{\psi^{N}}$

$\displaystyle\left\|{\Delta^{+}\ket{\psi^{N}}}\right\|^{2}\leq\left\|{\ket{\psi^{N}}}\right\|^{2}(N+1)d^{2}(v_{1},v_{2}),$
$\displaystyle\left\|{\Delta^{-}\ket{\psi^{N}}}\right\|^{2}\leq\left\|{\ket{\psi^{N}}}\right\|^{2}Nd^{2}(v_{1},v_{2}).$

Furthermore, the $N-$particle component of the state $\ket{\psi_{1}(t)}$ has a
bounded norm, with Eq. 8 providing the explicit bound. It thus immediately
follows that

$\displaystyle\left\|{\Delta^{+}\ket{\psi_{1}(t)}}\right\|^{2}\leq\big{(}1+\left\|{v_{1}}\right\|_{\infty}^{2}t\left\|{L}\right\|_{2}^{2}e^{2\pi\left\|{v_{1}}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}t}\big{)}d^{2}(v_{1},v_{2})$
$\displaystyle\left\|{\Delta^{-}\ket{\psi_{1}(t)}}\right\|^{2}\leq\left\|{v_{1}}\right\|_{\infty}^{2}t\left\|{L}\right\|_{2}^{2}e^{2\pi\left\|{v_{1}}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}t}d^{2}(v_{1},v_{2})$

These estimates, along with Eq. III.2, proves the statement of the lemma.
$\square$.

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_{\omega}\to
a_{\omega}e^{i\theta(\omega)}$ and $v(\omega)\to
v(\omega)e^{-i\theta(\omega)}$, for any $\theta:\mathbb{R}\to\mathbb{R}$,
leaves both the Hamiltonian in Eq. II.2 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.

###### Lemma 5

Let
$v_{1},v_{2}\in\textnormal{C}^{\infty}_{b}(\mathbb{R})\cap\textnormal{L}^{2}(\mathbb{R})$
be coupling functions and let $H_{1}(t),H_{2}(t)$ denote the Hamiltonians
corresponding to these coupling functions as given by Eq. II.2. Denoting by
$\rho_{i}(t)$ as the reduced state of the local system when an initial state
$\ket{\sigma}\otimes\ket{\textnormal{vac}}$, with
$\ket{\sigma}\in\mathcal{H}_{S}$, is evolved under the Hamiltonian $H_{i}(t),\
i\in\\{1,2\\}$, then

$\displaystyle\left\|{\rho_{1}(t)-\rho_{2}(t)}\right\|_{\textnormal{tr}}\leq\sqrt{2}g(t)\bigg{[}\int_{-\omega_{c}}^{\omega_{c}}\big{|}|v_{1}(\omega)|^{2}-|v_{2}(\omega)|^{2}\big{|}d\omega\bigg{]}^{1/2},$

where $g(t)$ is the function defined in lemma 4.

Proof: Denote by $\theta_{i}(\omega)$ the phase of the coupling coefficient
$v_{i}(\omega)\ \forall\omega\in\mathbb{R},i\in\\{1,2\\}$. Applying the
transformation $a_{\omega}\to a_{\omega}e^{-i\theta_{i}(\omega)}$ to the
Hamiltonian $H_{i}$, we obtain

$\displaystyle\tilde{H}_{i}(t)=$ $\displaystyle
H_{S}(t)+\int_{-\infty}^{\infty}\omega a_{\omega}^{\dagger}a_{\omega}d\omega+$
$\displaystyle\int_{-\infty}^{\infty}\big{|}v_{i}(\omega)\big{|}\big{(}a_{\omega}L^{\dagger}+a_{\omega}^{\dagger}L\big{)}d\omega.$

For $i\in\\{1,2\\}$, denote by $\ket{\psi_{i}(t)}$ the state obtained by
evolving $\ket{\sigma}\otimes\ket{\text{vac}}$ under the Hamiltonian
$\tilde{H}_{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\in\\{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
$\ket{{\psi_{1}(t)}}$ and $\ket{{\psi_{2}(t)}}$, which can be done with an
application of lemma 4:

$\displaystyle\left\|{\ket{\psi_{1}(t)}-\ket{\psi_{2}(t)}}\right\|\leq
g(t)d(|v_{1}|,|v_{2}|),$

where $g(t)$ and $d(v_{1},v_{2})$ are defined in lemma 4. Furthermore, we note
that

$d^{2}(|v_{1}|,|v_{2}|)\leq\int_{-\infty}^{\infty}\big{|}|v_{1}(\omega)|^{2}-|v_{2}(\omega)|^{2}\big{|}d\omega,$

and thus we obtain

$\left\|{\ket{\psi_{1}(t)}-\ket{\psi_{2}(t)}}\right\|\leq
g(t)\bigg{[}\int_{-\infty}^{\infty}\big{|}|v_{1}(\omega)|^{2}-|v_{2}(\omega)|^{2}\big{|}d\omega\bigg{]}^{1/2}$

Finally, using Eq. 3 together with the fact that partial trace is contractive,
we obtain the bound provided in the lemma $\square$.

We next consider the pseudomode approximation for non-Markovian dynamics the
coupling constant $v(\omega)$ by a sum of lorentzians, and solving a Markovian
master equation for the resulting approximation. This procedure, can be made
precise with the following result from Ref. Tamascelli  _et al._ (2018).

###### Lemma 6

[Pseudomodes from Ref. Tamascelli  _et al._ (2018)] Given a coupling constant
$|v(\omega)|^{2}$ which satisfies

$\displaystyle|v(\omega)|^{2}=\sum_{i=1}^{M}\frac{V_{i}}{(\omega-\omega_{i})^{2}+\kappa_{i}^{2}/4}\
\forall\ \omega\in\mathbb{R},$ (15)

where $V_{i}\geq 0$ and $\omega_{i}\in\mathbb{R}\ \forall i\in\\{1,2\dots
M\\}$, then the reduced state $\rho(t)$ of the local system obtained in
evolving an initial state $\ket{\sigma}\otimes\ket{\textnormal{vac}}$ under
the Hamiltonian given by Eq. II.2 can be computed exactly by solving the
following problem: Consider a system with Hilbert space
$\mathcal{H}_{S}\otimes\mathcal{H}_{\textnormal{aux}}$, where
$\mathcal{H}_{\textnormal{aux}}$ is the Fock space over $\mathbb{C}^{M}$, and
Hamiltonian

$\hat{H}(t)=H_{S}(t)+\sum_{i=1}^{M}\omega_{i}a_{i}^{\dagger}a_{i}+\sum_{i=1}^{M}\bigg{(}\frac{2\pi
V_{i}}{\kappa_{i}}\bigg{)}^{1/2}\big{(}a_{i}L^{\dagger}+La_{i}^{\dagger}\big{)},$

with $H_{S}(t)$ being the Hamiltonian of the local system introduced in Eq.
II.2 and $a_{1},a_{2}\dots a_{M}$ being the annihilation operators for the
Fock space $\mathcal{H}_{\textnormal{aux}}$. Compute the state $R(t)$ obtained
on evolving an initial state $\ket{\sigma}\otimes\ket{0}^{M}$ under the master
equation

$\displaystyle\frac{dR(t)}{dt}=-i[\hat{H}(t),R(t)]$
$\displaystyle\qquad+\kappa_{i}\sum_{i=1}^{M}\bigg{(}\frac{1}{2}\\{a_{i}^{\dagger}a_{i},R(t)\\}+a_{i}R(t)a_{i}^{\dagger}\bigg{)},$

then $\rho(t)=\textnormal{Tr}_{\textnormal{aux}}[R(t)]$.

Consequently, 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(\omega)|^{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
$\omega_{c}>0$, an integer $M$ indicating the number of Lorentzians and solve
the following optimization problem:

$\displaystyle\min_{\begin{subarray}{c}\omega\in\mathbb{R}^{M}\\\
V,\kappa\in\mathbb{R}^{+M}\end{subarray}}\int_{-\infty}^{\infty}\bigg{|}|v(\omega)|^{2}\textnormal{rect}_{\omega_{c}}(\omega)-\sum_{i=1}^{M}\frac{V_{i}}{(\omega-\omega_{i})^{2}+\kappa_{i}^{2}/4}\bigg{|}d\omega.$

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 (proved in appendix D), an upper bound
on the approximation error can be provided by an explicit Lorentzian
construction.

###### Lemma 7

Given a differentiable function $\Gamma:\mathbb{R}\to[0,\infty)$ which is also
bounded, then for $\omega_{c},\gamma\in(0,\infty)$, $M\in\mathbb{N}$,

$\displaystyle\int_{-\infty}^{\infty}\bigg{|}\Gamma(\omega)\textnormal{rect}_{\omega_{c}}(\omega)-\hat{\Gamma}(\omega)\bigg{|}d\omega\leq\varepsilon(\omega_{c},\kappa,M)$

where $\hat{\Gamma}:\mathbb{R}\to\mathbb{R}$ is a lorentzian approximation to
$\Gamma$ given by

$\hat{\Gamma}(\omega)=\sum_{n=1}^{M}\frac{\kappa\delta}{\pi}\frac{\Gamma(\omega_{n})}{(\omega-\omega_{n})^{2}+\kappa^{2}},$

with $\delta=2\omega_{c}/(M-1)$ and $\omega_{n}=-\omega_{c}+(n-1)\delta$ and
$\varepsilon(\omega_{c},\kappa,M)$ is the approximation error which satisfies
the scaling

$\displaystyle\varepsilon(\omega_{c},\kappa,M)=$
$\displaystyle\qquad\mathcal{O}\bigg{(}\Gamma_{\textnormal{max}}\kappa\log\frac{\kappa}{\omega_{c}}\bigg{)}+\mathcal{O}\bigg{(}\omega_{c}\kappa\Gamma^{d}_{\textnormal{max}}(\omega_{c})\log\frac{\kappa}{\omega_{c}}\bigg{)}+$
$\displaystyle\qquad\mathcal{O}\bigg{(}\frac{\omega_{c}^{3}}{\kappa
N}\Gamma^{d}_{\textnormal{max}}(\omega_{c})\bigg{)}+\mathcal{O}\bigg{(}\frac{\omega_{c}^{3}}{\kappa^{2}N}\Gamma_{\textnormal{max}}\bigg{)}+\mathcal{O}\bigg{(}\frac{\Gamma_{\textnormal{max}}}{M}\bigg{)},$

where $\Gamma_{\textnormal{max}}=\sup_{\omega\in\mathbb{R}}\Gamma(\omega)$ and
$\Gamma_{\textnormal{max}}^{d}(\omega_{c})=\sup_{\omega\in[-\omega_{c},\omega_{c}]}|\Gamma^{\prime}(\omega)|$.

Equipped with this error bound on the lorentzian approximation of the
magnitude square of the coupling function $v$, we can now show that the
pseudomode approximation is guaranteed to converge as well as provide a
scaling of the resulting approximation error with the number of pseudomodes.

Proof of theorem 2: For a coupling function
$v\in\text{C}^{\infty}_{b}(\mathbb{R})$, we denote by $\rho(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. II.2. First we introduce a
frequency cutoff $\omega_{c}$ — we denote by $\rho_{\omega_{c}}(t)$ the
reduced state of the local system obtained when evolved as per Hamiltonian in
Eq. 1. Second, we approximate $v\cdot\text{rect}_{\omega_{c}}$ by a sum of
lorentzians by using the construction of lemma 7 to obtain a pseudomode
approximation to the non-Markovian system (lemma 6). Let $\hat{\rho}(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
$\rho(t)$ and the reduced state using the pseudomode approximation
$\hat{\rho}(t)$ via

$\left\|{\rho(t)-\hat{\rho}(t)}\right\|_{\text{tr}}\leq\underbrace{\left\|{\rho(t)-\rho_{\omega_{c}}(t)}\right\|_{\text{tr}}}_{\text{Cutoff
error}}+\underbrace{\left\|{\rho_{\omega_{c}}(t)-\hat{\rho}(t)}\right\|_{\text{tr}}}_{\text{Lorentzian
approximation error}}.$

From theorem 1, it follows that the cutoff error decreases as
$\omega_{c}\to\infty$. To estimate the Lorentzian approximation error, we use
lemmas 5 and 7, which allows us to relate this error to the frequency cutoff
$\omega_{c}$, number of pseudomodes $M$ as well as the linewidths of the
pseudomodes $\kappa$. The key idea behind showing that the sum total of these
two errors goes to 0 as $M\to\infty$ is to choose $\kappa$ and $\omega_{c}$ as
a function of $M$ such that (i) the Lorentzian approximation error decreases
with $M$ and (ii) $\omega_{c}$ increases with $M$ (as a consequence of which
the cutoff error decreases with $M$). Below, we show that this choice can
always be made for any coupling function $v$.

We make the choice $\kappa=\Theta(M^{-1/4})$. Furthermore, consider the
increasing functions $f_{1},f_{2},f_{3},f_{4},f_{5}:[0,\infty)\to[0,\infty)$
of the cutoff frequencies that appear in the error estimate —
$f_{1}(\omega_{c})=\log\omega_{c}$,
$f_{2}(\omega_{c})=\omega_{c}\Gamma_{\text{max}}^{d}(\omega_{c})$,
$f_{3}=\omega_{c}\log(\omega_{c})\Gamma_{\text{max}}^{d}(\omega_{c})$,
$f_{4}(\omega_{c})=\omega_{c}^{3}$ and
$f_{5}(\omega_{c})=\omega_{c}^{3}\Gamma_{\text{max}}^{d}(\omega_{c})$. We note
that since by definition $\Gamma_{\text{max}}^{d}(\omega_{c})$ is an
increasing non-negative function, $f_{1},f_{2},f_{3},f_{4}=\mathcal{O}(f_{5})$
as $\omega_{c}\to\infty$. We choose $\omega_{c}=\Theta(f_{5}^{-1}(M^{1/5}))$ —
we note that since $f_{5}^{-1}$ is an increasing function,
$f_{1}(\omega_{c}),f_{2}(\omega_{c}),f_{3}(\omega_{c}),f_{4}(\omega_{c})=\mathcal{O}(M^{1/8})$
as $M\to\infty$. Using these choices of $\kappa$ and $\omega_{c}$, the error
in approximating $|v|^{2}$ with a sum of lorentzian, as calculated in Lemma 7,
is $\mathcal{O}(M^{-1/8})$ and thus goes to 0 as $M\to\infty$. Furthermore,
since $\omega_{c}$ grows with $M$, it follows from theorem 1 that the cutoff
error goes to 0 as $M\to\infty$.

To show that, under the assumption of a polynomial growth in
$|v^{\prime}(\omega)|$ with $\omega$ and a polynomial fall off of the cutoff
error with $\omega$, the total approximation error decreases polynomially with
the number of pseudomodes, we simply have to note that in the above proof,
$f_{5}^{-1}(M^{1/5})$ will be $\text{poly}(M)$ and consequently the
construction of theorem 3 yields a cutoff frequency $\omega_{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 $\omega_{c}$ the cutoff error also reduces polynomially
with the number of pseudomodes $\square$.

###  III.3 Convergence of Wilson’s star-to-delta transformation

Finally, we consider the Wilson’s star-to-delta transformation for simulating
non-Markovian quantum systems. This transformation maps the simulation of the
non-Markovian quantum system to a Hamiltonian simulation problem via the
Lanczos iteration.

We first introduce a frequency cutoff $\omega_{c}$ into the problem — in order
to apply Lanczos iteration, we need to define a Hamiltonian which would
generate the Lanczos chain and an initial mode which serves as the starting
point for the Lanczos iteration. The Hamiltonian generating the Lanczos chain
is simply the Hamiltonian of the environment with a frequency cutoff i.e.

$H_{\omega_{c}}^{E}=\int_{-\omega_{c}}^{\omega_{c}}\omega
a_{\omega}^{\dagger}a_{\omega}d\omega,$

and the initial mode $A_{0}$ is constructed from the coupling function $v$:

$A_{0}=\frac{\int_{-\omega_{c}}^{\omega_{c}}v(\omega)a_{\omega}d\omega}{\sqrt{\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}d\omega}}.$

The Lanczos iteration then generates a sequence of modes,
$A_{0},A_{1},A_{2},A_{3}\dots$ which can be expressed as

$A_{i}=\frac{\int_{-\omega_{c}}^{\omega_{c}}p_{i}(\omega)v(\omega)a_{\omega}d\omega}{\sqrt{\int_{-\omega_{c}}^{\omega_{c}}p_{i}^{2}(\omega)|v(\omega)|^{2}d\omega}},$

where $p_{i}(\omega)$ is the $i^{\text{th}}$ orthogonal polynomial generated
from $p_{0}(\omega)=1$ and which satisfy the recurrence

$\displaystyle
p_{i+1}(\omega)=(\omega-\alpha_{i})p_{i}(\omega)-\beta_{i}p_{i-1}(\omega),$
(16a) with $\beta_{0}=0$ and
$\displaystyle\alpha_{i}=\frac{\int_{-\omega_{c}}^{\omega_{c}}\omega
p_{i}^{2}(\omega)|v(\omega)|^{2}d\omega}{\int_{-\omega_{c}}^{\omega_{c}}p_{i}^{2}(\omega)|v(\omega)|^{2}d\omega}\
\text{and}\
\beta_{i}=\frac{\int_{-\omega_{c}}^{\omega_{c}}p_{i-1}^{2}(\omega)|v(\omega)|^{2}d\omega}{\int_{-\omega_{c}}^{\omega_{c}}p_{i}^{2}(\omega)|v(\omega)|^{2}d\omega}.$
(16b)

The operators $A_{i}$ satisfy $[A_{i},A_{j}^{\dagger}]=\delta_{i,j}$ and the
Hamiltonian $H_{\omega_{c}}^{E}$ can be expressed as

$\displaystyle
H_{\omega_{c}}^{E}=\sum_{i=0}^{\infty}\alpha_{i}A_{i}^{\dagger}A_{i}+\sum_{i=0}^{\infty}\sqrt{\beta_{i+1}}\big{(}A_{i}A_{i+1}^{\dagger}+A_{i+1}A_{i}^{\dagger}\big{)}.$
(17)

Furthermore, with this representation of the environment, the full local
system-environment Hamiltonian, with a frequency cutoff, can be expressed as

$\displaystyle
H_{\omega_{c}}(t)=H_{S}(t)+H_{\omega_{c}}^{E}+(L^{\dagger}B+B^{\dagger}L),$
(18)

where
$B=\int_{-\omega_{c}}^{\omega_{c}}v(\omega)a_{\omega}d\omega=\big{(}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}d\omega\big{)}^{1/2}A_{0}$.

The representation of the Hamiltonian in Eqs. 17 and 18 with the (infinite
sequence of) Lanczos modes is exact within the frequency cutoff. To perform a
numerical simulation, we need to approximate the infinite sequence of Lanczos
modes with a finite number of modes. With $N$ modes, the environment
Hamiltonian can be approximated as

$\displaystyle
H_{\omega_{c},N}^{E}=\sum_{i=0}^{N-1}\alpha_{i}A_{i}^{\dagger}A_{i}+\sum_{i=0}^{N-2}\sqrt{\beta_{i+1}}(A_{i}A_{i+1}^{\dagger}+A_{i+1}A_{i}^{\dagger}),$
(19)

and the full Hamiltonian can be approximated by

$\displaystyle
H_{\omega_{c},N}(t)=H_{S}(t)+H_{\omega_{c},N}^{E}+(L^{\dagger}B+B^{\dagger}L).$
(20)

In order to assess the convergence and accuracy of the approximation, it is of
interest to compare the quantum state obtained on separately evolving the
local system with $H_{\omega_{c}}(t)$ and $H_{\omega_{c},N}(t)$. It is easier
to do this in the interaction picture with respect to the environment
Hamiltonian, as made concrete in the following lemma.

###### Lemma 8

Let $\ket{\psi_{\omega_{c}}(t)}$ and $\ket{\psi_{\omega_{c},N}(t)}$ be the
joint state of the local system and environment obtained on evolving an
initial state $\ket{\sigma}\otimes\ket{\textnormal{vac}}$, with
$\ket{\sigma}\in\mathcal{H}_{S}$, with respect to $H_{\omega_{c}}(t)$ and
$H_{\omega_{c},N}(t)$ respectively, then

$\displaystyle\left\|{\rho_{\omega_{c}}(t)-{\rho_{\omega_{c},N}(t)}}\right\|_{\textnormal{tr}}$
$\displaystyle\qquad\leq\sqrt{2}g(t)\bigg{(}\sup_{s\in[0,t]}\Big{|}[\delta
B_{\omega_{c},N}(s),\delta
B_{\omega_{c},N}^{\dagger}(s)]\Big{|}\bigg{)}^{1/2},$

where

$\displaystyle
g(t)=\int_{0}^{t}\bigg{(}\sqrt{1+\left\|{v}\right\|_{\infty}^{2}\tau\left\|{L}\right\|_{2}^{2}e^{2\pi\left\|{v}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}\tau}}+$
$\displaystyle\left\|{v}\right\|_{\infty}\sqrt{\tau}\left\|{L}\right\|_{2}e^{\pi\left\|{v}\right\|^{2}_{\infty}\left\|{L}\right\|_{2}^{2}\tau}\bigg{)}d\tau,$

and

$\delta
B_{\omega_{c},N}(s)=e^{isH_{\omega_{c}}^{E}}Be^{-isH_{\omega_{c}}^{E}}-e^{isH_{\omega_{c},N}^{E}}Be^{-isH_{\omega_{c},N}^{E}}.$

Proof: The proof of this lemma is very similar to that of lemma 4, but
repeating the analysis in the interaction picture with respect to the
Hamiltonian of the environment ($H_{\omega_{c}}^{E}$ or $H_{\omega_{c},N}^{E}$
for the two models under consideration). The two Hamiltonians under
consideration, after moving into the interaction picture, denoted by
$\bar{H}_{\omega_{c}}(t)$ and $\bar{H}_{\omega_{c},N}(t)$ are given by

$\displaystyle\bar{H}_{\omega_{c}}(t)=H_{S}(t)+\big{(}B_{\omega_{c}}(t)L^{\dagger}+B^{\dagger}_{\omega_{c}}(t)L\big{)},$
$\displaystyle\bar{H}_{\omega_{c},N}(t)=H_{S}(t)+\big{(}B_{\omega_{c},N}(t)L^{\dagger}+B^{\dagger}_{\omega_{c},N}(t)L\big{)},$

where $B_{\omega_{c}}(t)=e^{itH_{\omega_{c}}^{E}}Be^{-itH_{\omega_{c}}}$ and
$B_{\omega_{c},N}(t)=e^{itH_{\omega_{c},N}^{E}}Be^{-itH_{\omega_{c},N}^{E}}$.
We denote by $\ket{\psi_{\omega_{c}}(t)}$ and $\ket{\psi_{\omega_{c},N}(t)}$
as the solution of the Schroedinger equation with Hamiltonians
$\tilde{H}_{\omega_{c}}(t)$ and $\tilde{H}_{\omega_{c},N}(t)$ and with the
initial state being $\ket{\sigma}\otimes\ket{\text{vac}}$. We note that the
norm of the $N-$particle component of $\ket{\psi_{\omega_{c}}(t)}$ will
satisfy the bound in Eq. 8 since it corresponds to the quantum state obtained
from a Hamiltonian of the form Eq. II.2 with coupling function
$v\cdot\text{rect}_{\omega_{c}}$, and translating this state into the
interaction picture does not change these norms. Consequently, we obtain

$\displaystyle\left\|{\ket{\psi_{\omega_{c},N}(t)}-\ket{\psi_{\omega_{c}}(t)}}\right\|$
$\displaystyle\leq\int_{0}^{t}\left\|{(\bar{H}_{\omega_{c}}(s)-\bar{H}_{\omega_{c},N}(s))\ket{\psi_{\omega_{c}}(s)}}\right\|ds,$
$\displaystyle\leq\left\|{L}\right\|_{2}\int_{0}^{t}\bigg{(}\left\|{\delta
B_{\omega_{c},N}(s)\ket{\psi_{\omega_{c}}(s)}}\right\|+\left\|{\delta
B_{\omega_{c},N}^{\dagger}(s)\ket{\psi_{\omega_{c}}(s)}}\right\|\bigg{)}ds.$

Noting that for an $N-$particle wave-packet $\ket{\psi^{N}}$ and for all
$s\in[0,t]$

$\displaystyle\left\|{\delta B_{\omega_{c},N}(s)\ket{\psi^{N}}}\right\|^{2}$
$\displaystyle\qquad\leq\left\|{\ket{\psi^{N}}}\right\|^{2}N\sup_{s\in[0,t]}\Big{|}[\delta
B_{\omega_{c},N}(s),\delta B_{\omega_{c},N}^{\dagger}(s)]\Big{|},$
$\displaystyle\left\|{\delta
B_{\omega_{c},N}^{\dagger}(s)\ket{\psi^{N}}}\right\|^{2}$
$\displaystyle\qquad\leq\left\|{\ket{\psi^{N}}}\right\|^{2}(N+1)\sup_{s\in[0,t]}\Big{|}[\delta
B_{\omega_{c},N}(s),\delta B_{\omega_{c},N}^{\dagger}(s)]\Big{|}.$

Using this together with the bound in Eq. 8, we obtain that

$\displaystyle\left\|{\ket{\psi_{\omega_{c}}(t)}-\ket{\psi_{\omega_{c},N}(t)}}\right\|$
$\displaystyle\qquad\leq g(t)\bigg{(}\sup_{s\in[0,t]}\Big{|}[\delta
B_{\omega_{c},N}(s),\delta
B_{\omega_{c},N}^{\dagger}(s)]\Big{|}\bigg{)}^{1/2}.$

Finally, the lemma follows using the contractivity of the partial trace
together with Eq. 3 $\square$.

An explicit dependence of the error bound in lemma 8 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 and proved explicitly in appendix E.

###### Lemma 9

Let $\delta B_{\omega_{c},N}(s)$ be as defined in lemma 8, then

$\Big{|}[\delta B_{\omega_{c},N}(s),\delta
B_{\omega_{c},N}^{\dagger}(s)]\Big{|}\leq\mathcal{O}\bigg{(}\frac{\left\|{v}\right\|_{\infty}^{2}(2\omega_{c})^{N+2}s^{N+1}}{(N-1)!}\bigg{)}.$

Equipped with this error bound, we can now show that for sufficiently large
frequency cutoff and sufficiently large number of modes, the Wilson star-to-
delta transformation is guaranteed to converge to the correct result.

Proof of theorem 3: The proof of this theorem parallels that of theorem 2. For
a coupling function $v\in\text{C}^{\infty}_{b}(\mathbb{R})$, we denote by
$\rho(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. II.2. First
we introduce a frequency cutoff $\omega_{c}$ — we denote by
$\rho_{\omega_{c}}(t)$ the reduced state of the local system obtained when
evolved as per Hamiltonian in Eq. 1. Second, we use the Wilson’s star-to-delta
transformation with $M$ modes to obtain an approximation $\hat{\rho}(t)$ to
$\rho_{\omega_{c}}(t)$. We estimate the error between the true reduced state
$\rho(t)$ and the reduced state using the pseudomode approximation
$\hat{\rho}(t)$ via

$\left\|{\rho(t)-\hat{\rho}(t)}\right\|_{\text{tr}}\leq\underbrace{\left\|{\rho(t)-\rho_{\omega_{c}}(t)}\right\|_{\text{tr}}}_{\text{Cutoff
error}}+\underbrace{\left\|{\rho_{\omega_{c}}(t)-\hat{\rho}(t)}\right\|_{\text{tr}}}_{\text{Wilson
star-to-delta error}}.$

From theorem 1, it follows that the cutoff error decreases as
$\omega_{c}\to\infty$. To estimate
$\left\|{\rho_{\omega_{c}}(t)-\hat{\rho}(t)}\right\|_{\text{tr}}$, we use
lemmas 8 and 9, which allows us to relate this error to the frequency cutoff
$\omega_{c}$ and the number of modes $M$. We note that choosing $\omega_{c}$
to increase sublinearly with $M$, both the cutoff error
$\left\|{\rho(t)-\rho_{\omega_{c}}(t)}\right\|_{\text{tr}}$ and
$\left\|{\rho_{\omega_{c}}(t)-\hat{\rho}(t)}\right\|_{\text{tr}}$ go to 0 as
$M\to\infty$, thus providing a convergence guarantee for Wilson star-to-delta
transformation. Furthermore, with this choice of $\omega_{c}$ and under the
assumption that the cutoff error falls polynomially with $\omega_{c}$, it
immediately follows that the approximation error decreases polynomially with
the number of modes $M$ $\square$.

##  IV Conclusions and outlook

Our work provides a rigorous analysis of the pseudomode approximation to
simulating 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. Our results
provide a rigorous footing to numerical methods for simulating non-Markovian
quantum systems used in practice.

Several questions broadly 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 for which the resulting system-environment
dynamics is well defined. For instance, a number of models in practice use a
point-coupling model in between the local system and a gaussian bosonic
environment which cannot be described by a square integrable coupling function
yet generate physically reasonable dynamics. The theory of single point
coupling (which generates Markovian local system dynamics) can be rigorously
dealt with using the machinery of quantum stochastic calculus, however it is
unclear if extensions to multiple point couplings (which generates non-
Markovian local system dynamics) is possible. 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
made in section II.2.

Second, as we mentioned in the previous sections, our error estimates grow
exponentially with time. This appears to be a consequence of the fact that the
local system can emit arbitrarily large number of excitations into the
environment, and the only generally provable bound on the probability of
emission of $N-$particles is poissonian in $N$ with parameter that grows with
$t$. However, for practically useful problems, we expect that the
approximation error will grow only polynomially with time (for e.g. if the
environment was finite dimensional, it can be trivially shown that the error
grows linearly with time) and leave a possible improvement in this bound for
future work.

###### Acknowledgements.

RT acknowledges Max Planck Harvard research center for quantum optics (MPHQ)
postdoctoral fellowship. We acknowledge support from the ERC Advanced Grant
QUENOCOBA under the EU Horizon 2020 program (grant agreement 742102) and from
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
the project number414325145 in the framework of the Austrian Science
Fund(FWF): SFB F7104.

##  Appendix A Systems for which assumption 1 is provable

###### Proposition 1

Assumption 1 is true for coupling functions $v$ that are bounded and square
integrable,
$v\in\textnormal{B}(\mathbb{R})\cap\textnormal{L}^{2}(\mathbb{R})$.

Proof: _Square-integrable coupling constants:_ We consider constants
$v:\mathbb{R}\to\mathbb{C}$ that are smooth, bounded and square integrable. We
note that in this case, for a symmetric and positive kernel
$K:[0,t]^{2}\to\mathbb{C}$,

$\displaystyle\lim_{\Omega\to\infty}\int_{\omega_{c}\leq|\omega_{1}|,|\omega_{2}|\leq\Omega}|v(\omega_{1})|^{2}|v(\omega_{2})|^{2}\bigg{(}\int_{s_{1},s_{2}=0}^{t}K(s_{1},s_{2})e^{i\omega_{2}(t-s_{2})}e^{-i\omega_{1}(t-s_{1})}ds\bigg{)}d\omega_{1}d\omega_{2}=$
$\displaystyle\bigg{|}\int_{\omega_{c}\leq|\omega_{1}|,|\omega_{2}|}|v(\omega_{1})|^{2}|v(\omega_{2})|^{2}\bigg{(}\int_{s_{1},s_{2}=0}^{t}K(s_{1},s_{2})e^{i\omega_{2}(t-s_{2})}e^{-i\omega_{1}(t-s_{1})}ds\bigg{)}d\omega_{1}d\omega_{2}\bigg{|}\leq$
$\displaystyle\int_{\omega_{c}\leq|\omega_{1}|,|\omega_{2}|}|v(\omega_{1})|^{2}|v(\omega_{2})|^{2}\bigg{(}\int_{s_{1},s_{2}=0}^{t}|K(s_{1},s_{2})|ds_{1}ds_{2}\bigg{)}d\omega_{1}d\omega_{2}\leq
t^{2}\bigg{(}\int_{\omega_{c}\leq|\omega|}|v(\omega)|^{2}d\omega\bigg{)}^{2}\sup_{s_{1},s_{2}\in[0,t]}|K(s_{1},s_{2})|.$

By noting that
$\sup_{s_{1},s_{2}\in[0,t]}|K(s_{1},s_{2})|\leq\left\|{K}\right\|_{[0,t]^{2}}^{\mathcal{S}}$,
it follows that square-integrable coupling constants satisfy the conditions of
assumption 1 with $V(\omega_{c},t)$ being given by

$V(\omega_{c},t)=t^{2}\int_{|\omega|\geq\omega_{c}}|v(\omega)|^{2}d\omega\
\qquad\qquad\qquad\qquad\square$

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.

###### Lemma 10

Let $f\in\textnormal{AC}([0,a])$, then $\forall p>0$

$\bigg{|}\frac{\pi}{2}f(0)-\int_{0}^{a}\frac{\sin
px}{x}f(x)dx\bigg{|}\leq\mathcal{O}\bigg{(}\frac{1}{\sqrt{ap}}\bigg{)}\left\|{f}\right\|^{\mathcal{S}}_{[0,a]}.$

Proof: We begin by using the Dirchlet integral, $\int_{0}^{\infty}\sin
x/xdx=\pi/2$, from which it follows that

$\bigg{|}\frac{\pi}{2}f(0)-\int_{0}^{a}\frac{\sin
px}{x}f(x)dx\bigg{|}\leq\bigg{|}\int_{0}^{a}\big{(}f(x)-f(0)\big{)}\frac{\sin
px}{x}dx\bigg{|}+|f(0)|\bigg{|}\int_{a}^{\infty}\frac{\sin px}{x}dx\bigg{|}.$

The second term is easily bounded

$\displaystyle|f(0)|\bigg{|}\int_{a}^{\infty}\frac{\sin px}{x}dx\bigg{|}$
$\displaystyle\leq\bigg{(}\sup_{x\in[0,a]}|f(x)|\bigg{)}\bigg{|}-\frac{\cos
pa}{pa}+\int_{a}^{\infty}\frac{\cos px}{px^{2}}dx\bigg{|}$
$\displaystyle\leq\bigg{(}\sup_{x\in[0,a]}|f(x)|\bigg{)}\bigg{(}\frac{1}{pa}+\int_{a}^{\infty}\frac{1}{px^{2}}dx\bigg{)}=\bigg{(}\sup_{x\in[0,a]}|f(x)|\bigg{)}\frac{2}{pa}.$

To bound the first term, we split the integral at $b(p)\in[0,a]$, where
$b(p)$, to be specified later, vanishes as $p\to\infty$:

$\bigg{|}\int_{0}^{a}\big{(}f(x)-f(0)\big{)}\frac{\sin
px}{x}dx\bigg{|}\leq\int_{0}^{b(p)}\big{|}f(x)-f(0)\big{|}\frac{dx}{x}+\bigg{|}\int_{b(p)}^{a}\big{(}f(x)-f(0)\big{)}\frac{\sin
px}{x}dx\bigg{|}.$

Noting that $|f(x)-f(0)|\leq|x|\ \text{esssup}_{x\in[0,a]}|f^{\prime}(x)|$, we
obtain

$\int_{0}^{b(p)}\big{|}f(x)-f(0)\big{|}\frac{dx}{x}\leq b(p)\
\underset{x\in[0,a]}{\text{ esssup }}|f^{\prime}(x)|.$

Futhermore, using integration by parts, we obtain

$\displaystyle\bigg{|}\int_{b(p)}^{a}\big{(}f(x)-f(0)\big{)}\frac{\sin
px}{x}dx\bigg{|}$
$\displaystyle=\bigg{|}\frac{f(b(p))\cos(pb(p))}{pb(p)}-\frac{f(a)\cos
pa}{pa}+\frac{1}{p}\int_{b(p)}^{a}\bigg{(}\frac{f^{\prime}(x)}{x}-\frac{f(x)-f(0)}{x^{2}}\bigg{)}\cos
pxdx\bigg{|},$
$\displaystyle\leq\frac{3}{p}\bigg{(}\frac{1}{b(p)}+\frac{1}{a}\bigg{)}\sup_{x\in[0,a]}|f(x)|+\frac{\log(a/b(p))}{p}\underset{x\in[0,a]}{\text{
esssup }}|f^{\prime}(x)|.$

Using these estimates and noting that
$\sup_{x\in[0,a]}|f(x)|\leq\left\|{f}\right\|_{[0,a]}^{\mathcal{S}}$ and
$\text{esssup}_{x\in[0,a]}|f^{\prime}(x)|\leq\left\|{f}\right\|_{[0,a]}^{\mathcal{S}}/a$,
we obtain

$\bigg{|}\frac{\pi}{2}f(0)-\int_{0}^{a}\frac{\sin
px}{x}f(x)dx\bigg{|}\leq\bigg{(}\frac{b(p)}{a}+\frac{1}{p}\bigg{(}\frac{3}{b(p)}+\frac{5}{a}\bigg{)}+\frac{\log(a/b(p))}{ap}\bigg{)}\left\|{f}\right\|_{[0,a]}^{\mathcal{S}}.$

Since this bound holds for any $b(p)\in[0,a]$, we choose
$b(p)=\Theta(\sqrt{a/p})$ which proves the lemma $\square$.

###### Lemma 11

Let $f\in\textnormal{AC}([a,b])$ for $0<a<b$, then $\forall p>0$

$\bigg{|}\int_{a}^{b}\frac{\sin
px}{x}f(x)dx\bigg{|}\leq\mathcal{O}\bigg{(}\frac{\log(b/a)}{ap}\bigg{)}\left\|{f}\right\|_{[a,b]}^{\mathcal{S}}$

Proof: This follows straightforwardly from integration by parts.

$\int_{a}^{b}\frac{\sin px}{x}f(x)dx=\frac{f(a)\cos pa}{pa}-\frac{f(b)\cos
pb}{pb}+\frac{1}{p}\int_{a}^{b}\bigg{(}\frac{f^{\prime}(x)}{x}-\frac{f(x)}{x^{2}}\bigg{)}\cos
px\ dx,$

and consequently

$\bigg{|}\int_{a}^{b}\frac{\sin
px}{x}f(x)dx\bigg{|}\leq\bigg{(}\frac{2}{pa}+\frac{2}{pb}+\frac{\log(b/a)}{ap}\bigg{)}\left\|{f}\right\|^{\mathcal{S}}_{[a,b]},$

from which the lemma statement follows $\square$.

###### Lemma 12

Let $f\in\textnormal{AC}([a,b])$, then for $p>0$ and $y\in\mathbb{R}$

$\bigg{|}\int_{a}^{b}\bigg{(}\pi\delta(x-y)-\frac{\sin
p(x-y)}{x-y}\bigg{)}f(x)dx\bigg{|}\leq\xi(a,b,y,p)\left\|{f}\right\|_{[a,b]}^{\mathcal{S}}.$

where $\xi(a,b,y,p)\to 0$ as $p\to\infty$ and $\delta(\cdot)$ is the dirac-
delta function, whose action on a continuous function $f$ within an interval
$[a,b]$ is given by

$\int_{a}^{b}\delta(x-y)f(x)dx=\begin{cases}f(y)&\text{if }y\in(a,b)\\\
f(y)/2&\text{if }y\in\\{a,b\\}\\\ 0&\text{otherwise}\end{cases}.$

Proof: We consider three different cases

  1. 1.

If $y\in\\{a,b\\}$, then it follows from lemma 10 that

$\bigg{|}\int_{a}^{b}\bigg{(}\pi\delta(x-y)-\frac{\sin
p(x-y)}{x-y}\bigg{)}f(x)dx\bigg{|}\leq\mathcal{O}\bigg{(}\frac{1}{\sqrt{|a-b|p}}\bigg{)}\left\|{f}\right\|_{[a,b]}^{\mathcal{S}}.$

  2. 2.

If $y<a$ or $y>b$, then it follows from lemma 11 that

$\bigg{|}\int_{a}^{b}\bigg{(}\pi\delta(x-y)-\frac{\sin
p(x-y)}{x-y}\bigg{)}f(x)dx\bigg{|}\leq\mathcal{O}\bigg{(}\frac{|\log(|y-a|/|y-b|)|}{p\min(|y-a|,|y-b|)}\bigg{)}\left\|{f}\right\|_{[a,b]}^{\mathcal{S}}.$

  3. 3.

If $y\in(a,b)$, then

$\displaystyle\bigg{|}\int_{a}^{b}\bigg{(}\pi\delta(x-y)-\frac{\sin
p(x-y)}{x-y}\bigg{)}f(x)dx\bigg{|}$
$\displaystyle\leq\bigg{|}\frac{\pi}{2}f(y)-\int_{a}^{y}\frac{\sin
p(x-y)}{x-y}f(x)dx\bigg{|}+\bigg{|}\frac{\pi}{2}f(y)-\int_{y}^{b}\frac{\sin
p(x-y)}{x-y}f(x)dx\bigg{|}$
$\displaystyle\leq\mathcal{O}\bigg{(}\frac{1}{\sqrt{|y-a|p}}\bigg{)}\left\|{f}\right\|_{[a,y]}^{\mathcal{S}}+\mathcal{O}\bigg{(}\frac{1}{\sqrt{|y-b|p}}\bigg{)}\left\|{f}\right\|_{[y,b]}^{\mathcal{S}}$
$\displaystyle\leq\mathcal{O}\bigg{(}\frac{1}{\sqrt{\min(|y-a|,|y-b|)p}}\bigg{)}\left\|{f}\right\|_{[a,b]}^{\mathcal{S}}.$

Here we have used lemma 10 and the fact that $\forall y\in(a,b):\
\left\|{f}\right\|_{[a,y]}^{\mathcal{S}}\leq\left\|{f}\right\|_{[a,b]}^{\mathcal{S}}$
and
$\left\|{f}\right\|_{[y,b]}^{\mathcal{S}}\leq\left\|{f}\right\|_{]a,b]}^{\mathcal{S}}$
$\square$.

###### Lemma 13

For $K\in\textnormal{AC}_{\textnormal{sym}}([a,b]^{2})$ and
$y_{1},y_{2}\in\mathbb{R}$, then $\forall p>0$

$\displaystyle\bigg{|}\int_{x_{1},x_{2}=a}^{b}\bigg{(}\pi\delta(x_{1}-y_{1})-\frac{\sin
p(x_{1}-y_{1})}{x_{1}-y_{1}}\bigg{)}\bigg{(}\pi\delta(x_{2}-y_{2})-\frac{\sin
p(x_{2}-y_{2})}{x_{2}-y_{2}}\bigg{)}K(x_{1},x_{2})dx_{1}dx_{2}\bigg{|}$
$\displaystyle\leq\xi(a,b,y_{1},p)\xi(a,b,y_{2},p)\left\|{K}\right\|_{[a,b]^{2}}^{\mathcal{S}}.$

Proof: For $y\in[a,b]$, we define $\varphi_{y}:[a,b]\to\mathbb{C}$ via

$\varphi_{y}(x)=\int_{a}^{b}\bigg{(}\pi\delta(x^{\prime}-y)-\frac{\sin
p(x^{\prime}-y)}{x^{\prime}-y}\bigg{)}K(x,x^{\prime})dx^{\prime},$

such that

$\displaystyle\int_{x_{1},x_{2}=a}^{b}\bigg{(}\pi\delta(x_{1}-y_{1})-\frac{\sin
p(x_{1}-y_{1})}{x_{1}-y_{1}}\bigg{)}\bigg{(}\pi\delta(x_{2}-y_{2})-\frac{\sin
p(x_{2}-y_{2})}{x_{2}-y_{2}}\bigg{)}K(x_{1},x_{2})dx_{1}dx_{2}=$
$\displaystyle\int_{a}^{b}\bigg{(}\pi\delta(x_{1}-y_{1})-\frac{\sin
p(x_{1}-y_{1})}{x_{1}-y_{1}}\bigg{)}\varphi_{y_{2}}(x_{1})dx_{1}.$

Since $\varphi_{y}$ is absolutely continuous, it follows from lemma 12 that

$\displaystyle\bigg{|}\int_{a}^{b}\bigg{(}\pi\delta(x_{1}-y_{1})-\frac{\sin
p(x_{1}-y_{1})}{x_{1}-y_{1}}\bigg{)}\varphi_{y_{2}}(x_{1})dx_{1}\bigg{|}\leq\xi(a,b,y_{1},p)\left\|{\varphi_{y_{2}}}\right\|_{[a,b]}^{\mathcal{S}}.$
(21)

Furthermore, it also follows from lemma 12 and the fact that $K$ was
absolutely continuous in either of its arguments that

$\sup_{x\in[a,b]}|\varphi_{y_{2}}(x)|\leq\xi(a,b,y_{2},p)\left\|{K}\right\|_{[a,b]^{2}}^{\mathcal{S}}.$

Similarly, for $x$ almost everywhere in $[a,b]$, $\varphi_{y_{2}}^{\prime}(x)$
exists and consequently from lemma 12 it follows that

$\text{For }x\in_{a.e.}[a,b],\
(b-a)|\varphi^{\prime}_{y_{2}}(x)|\leq\xi(a,b,y_{2},p)\left\|{K}\right\|_{[a,b]^{2}}^{\mathcal{S}}.$

Consequently, $\left\|{\varphi_{y_{2}}}\right\|^{\mathcal{S}}_{[a,b]}\leq
2\xi(a,b,y_{2},p)\left\|{K}\right\|_{[a,b]^{2}}^{\mathcal{S}}$. Using this
along with Eq. 21 completes the proof $\square$.

###### Proposition 2

For $V_{1},V_{2}\dots V_{N}\in\mathbb{C}$,
$\tau_{1},\tau_{2}\dots\tau_{N}\in\mathbb{R}$, define a coupling function
$v\in\textnormal{C}_{b}^{\infty}(\mathbb{R})$ by
$v(\omega)=\sum_{i=1}^{N}V_{i}e^{i\omega\tau_{i}}\
\forall\omega\in\mathbb{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:

$|v(\omega)|^{2}=\sum_{i=1}^{\tilde{N}}\tilde{V}_{i}e^{i\omega T_{i}}.$

For any $K\in\text{AC}_{\text{sym}}([0,t]^{2})$, we obtain

$\displaystyle\lim_{\Omega\to\infty}\int_{\omega_{c}\leq|\omega_{1}|,|\omega_{2}|\leq\Omega}|v(\omega_{1})|^{2}|v(\omega_{2})|^{2}\bigg{(}\int_{s_{1},s_{2}=0}^{t}K(s_{1},s_{2})e^{i\omega_{2}(t-s_{2})}e^{-i\omega_{1}(t-s_{1})}ds\bigg{)}d\omega_{1}d\omega_{2}$
$\displaystyle=4\bigg{|}\sum_{i,j=1}^{\tilde{N}}\tilde{V}_{i}\tilde{V}_{j}^{*}\int_{s_{1},s_{2}=0}^{t}\bigg{(}\pi\delta(s_{2}-t-T_{i})-\frac{\sin\omega_{c}(s_{2}-t-T_{i})}{s_{2}-t-T_{i}}\bigg{)}\bigg{(}\pi\delta(s_{1}-t-T_{j})-\frac{\sin\omega_{c}(s_{1}-t-T_{j})}{s-t-
T_{j}}\bigg{)}K(s_{1},s_{2})ds_{1}ds_{2}\bigg{|}$ $\displaystyle\leq
4\sum_{i,j=1}^{\tilde{N}}|\tilde{V}_{i}\tilde{V}_{j}^{*}|\bigg{|}\int_{s_{1},s_{2}=0}^{t}\bigg{(}\pi\delta(s_{2}-t-T_{i})-\frac{\sin\omega_{c}(s_{2}-t-T_{i})}{s_{2}-t-T_{i}}\bigg{)}\bigg{(}\pi\delta(s_{1}-t-T_{j})-\frac{\sin\omega_{c}(s_{1}-t-T_{j})}{s-t-
T_{j}}\bigg{)}K(s_{1},s_{2})ds_{1}ds_{2}\bigg{|}.$

Applying lemma 13, we immediately obtain

$\displaystyle\lim_{\Omega\to\infty}\int_{\omega_{c}\leq|\omega_{1}|,|\omega_{2}|\leq\Omega}|v(\omega_{1})|^{2}|v(\omega_{2})|^{2}\bigg{(}\int_{s_{1},s_{2}=0}^{t}K(s_{1},s_{2})e^{i\omega_{2}(t-s_{2})}e^{-i\omega_{1}(t-s_{1})}ds\bigg{)}d\omega_{1}d\omega_{2}$
$\displaystyle\qquad\leq\bigg{(}\sum_{i,j=1}^{\tilde{N}}8\big{|}\tilde{V}_{i}\tilde{V}_{j}\big{|}\xi(0,t,t+T_{i},\omega_{c})\xi(0,t,t+T_{j},\omega_{c})\bigg{)}\left\|{K}\right\|_{[0,t]^{2}}^{\mathcal{S}},$

which proves that assumption 1 is satisfied $\square$.

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.

##  Appendix B $N$-point Green’s function and the environment state

In this appendix, we relate the $N-$point Green’s function defined in section
II.2 to the joint state of the local system and environment. We assume that
the initial state to be $\ket{\sigma}\otimes\ket{\text{vac}}$, where
$\ket{\sigma}\in\mathcal{H}_{S}$ and let
$\ket{\psi(t)}=U(t,0)\ket{\sigma}\otimes\ket{\text{vac}}$, where
$U(\tau_{1},\tau_{2})$ is the propagator corresponding to the Hamiltonian in
Eq. II.2. It immediately follows that

$\displaystyle\ket{\psi(t)}=\sum_{N=0}^{\infty}\ket{\psi^{N}(t)}=\sum_{N=0}^{\infty}\sum_{i=1}^{D}\frac{1}{N!}\int_{\omega\in\mathbb{R}^{N}}\psi^{N}_{i}(\omega;t)\bigg{(}\prod_{j=1}^{N}a_{\omega_{j}}^{\dagger}\bigg{)}\ket{e_{i}}\otimes\ket{\text{vac}},$
(22)

where $\\{\ket{e_{i}},i\in\\{1,2\dots D\\}\\}$ is a basis for the
$D-$dimensional local system and

$\psi^{N}_{i}(\omega;t)=\bra{e_{i}}\otimes\bra{\text{vac}}\bigg{(}\prod_{i=1}^{M}a_{\omega_{i}}\bigg{)}U(t,0)\ket{\sigma}\otimes\ket{\text{vac}}=\bra{e_{i}}\otimes\bra{\text{vac}}U(t,0)\mathcal{T}\bigg{(}\prod_{i=1}^{M}a_{\omega_{i}}(t)\bigg{)}\ket{\sigma}\otimes\ket{\text{vac}},$

where $a_{\omega}(t)=U(0,t)a_{\omega}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

$a_{\omega}(t)=a_{\omega}(0)e^{-i\omega
t}-iv^{*}(\omega)\int_{0}^{t}L(s)e^{-i\omega(t-s)}ds,$

from which it immediately follows that

$\displaystyle\psi^{N}_{i}(\omega;t)=(-i)^{N}\bigg{(}\prod_{j=1}^{N}v^{*}(\omega_{j})\bigg{)}\bra{e_{i}}\otimes\bra{\text{vac}}F_{N}(\omega;t)\ket{\sigma}\otimes\ket{\text{vac}},$
(23)

where $F_{N}(\omega;t)$ is defined in Eq. 7. Substituting Eq. 23 into Eq. 22
and using $\sum_{i=1}^{D}\ket{e_{i}}\bra{e_{i}}=\text{I}$, we obtain the
result in Eq. 7.

##  Appendix C Settings where assumption 2 is provable

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.

###### Proposition 3

Consider a coupling function $v$ which is a constant i.e.
$\forall\omega\in\mathbb{R}:v(\omega)=v_{0}$ for some constant $v_{0}$, then
assumption 2 is satisfied.

Proof: We begin by noting that constant coupling functions imply Markovian
dynamics for the local system, in which case the quantum regression theorem
implies that the $N-$point Green’s function can be expressed entirely evolving
the local system under a non-Hermitian effective Hamiltonian. More
specifically,

$G_{N}(s_{1},s_{2}\dots
s_{N};t)=U_{\text{eff}}(t,0)\mathcal{T}\big{(}L_{\text{eff}}(s_{1})L_{\text{eff}}(s_{2})\dots
L_{\text{eff}}(s_{N})\big{)},$

where $U_{\text{eff}}(\tau_{1},\tau_{2})$ is the propagator corresponding to
the Hamiltonian $H_{\text{eff}}(t)=H_{S}(t)-i\pi|v_{0}|^{2}L^{\dagger}L/2$ and
$L_{\text{eff}}(s)=U_{\text{eff}}(0,s)LU_{\text{eff}}(s,0)$. We note that
$\partial_{s}L_{\text{eff}}(s)=-i[L_{\text{eff}}(s),H_{\text{eff}}(s)]$ and
consequently $\forall s\in[0,t]:\
\left\|{\partial_{s}L_{\text{eff}}(s)}\right\|\leq\sup_{s\in[0,t]}\left\|{[L_{\text{eff}}(s),H_{\text{eff}}(s)]}\right\|_{2}$.
Noting that for $s_{2},s_{3}\dots s_{N}\in[0,t]$, as a function of
$s_{1}\in[0,t]$ $G_{N}(s_{1},s_{2}\dots s_{N};t)$ is differentiable at all
points except at $s_{1}\in\\{s_{2},s_{3}\dots s_{N}\\}$, we obtain

$\left\|{\partial_{s}G_{N}(s,s_{2}\dots
s_{N};t)}\right\|_{2}\leq\left\|{L}\right\|_{2}^{N-1}\sup_{s\in[0,t]}\left\|{[L_{\text{eff}}(s),H_{\text{eff}}(s)}\right\|_{2}\
\text{for }s\in_{a.e.}[0,t],$

wherein we have used that since
$\left\|{U_{\text{eff}}(\tau_{1},\tau_{2})}\right\|_{2}\leq 1$,
$\left\|{L_{\text{eff}}(s)}\right\|_{2}\leq\left\|{L}\right\|_{2}\ \forall\
s\in[0,t]$ $\square$.

###### Proposition 4

Consider a coupling function
$v\in\textnormal{C}_{b}^{\infty}(\mathbb{R})\cap\textnormal{L}^{2}(\mathbb{R})$,
then assumption 2 is satisfied.

Proof: Let $s\in[0,t]^{N-1}$ and without loss of generality, we assume that,
$s_{1}\leq s_{2}\dots\leq s_{N-1}$ since the $N-$point Green’s function is a
symmetric function of its arguments. We consider the $N-$point Green’s
function $G_{N}(\\{s_{0},s\\};t)$ as a function of $s_{0}$. Consider the case
where for some $i\in\\{0,1\dots N-2\\},s\in(s_{i},s_{i+1})$ (where
$s_{0}\equiv-\infty$ and $s_{N}\equiv\infty$) in which case the explicitly
applying the time-ordering operator we obtain

$G_{N}(\\{s_{0},s\\};t)=\bra{\text{vac}}U(t,0)L(s_{N-1})\dots
L(s_{i+1})L(s_{0})L(s_{i})\dots L(s_{1})\ket{\text{vac}}$

and therefore $\forall\ket{\sigma}\in\mathcal{H}_{S}$

$\left\|{\partial_{s_{0}}G_{N}(\\{s_{0},s\\};t)\ket{\sigma}}\right\|_{2}\leq\left\|{L}\right\|_{2}^{N-i-1}\left\|{[L,H(s_{0})]U(s_{0},0)L(s_{i})\dots
L(s_{1})\ket{\sigma}\otimes\ket{\text{vac}}}\right\|.$

Noting that

$[L,H(s_{0})]=[L,H_{S}(s_{0})]+[L,L^{\dagger}]\int_{-\infty}^{\infty}v(\omega)a_{\omega}d\omega,$

we obtain

$\displaystyle\left\|{\partial_{s_{0}}G_{N}(\\{s_{0},s\\};t)\ket{\sigma}}\right\|_{2}\leq$
$\displaystyle\qquad\left\|{L}\right\|_{2}^{N-1}\left\|{[L,H_{S}(s)]}\right\|_{2}+\left\|{[L,L^{\dagger}]}\right\|_{2}\left\|{L}\right\|_{2}^{N-i-1}\left\|{\bigg{(}\int_{-\infty}^{\infty}v(\omega)a_{\omega}d\omega\bigg{)}U(s_{0},0)L(s_{i})\dots
L(s_{1})\ket{\sigma}\otimes\ket{\text{vac}}}\right\|_{2}.$

To proceed further, we provide a bound on the second term in the above
estimate. Denote by $\ket{\phi}=U(s,0)L(s_{i})L(s_{i-1})\dots
L(s_{1})\ket{\sigma}\otimes\ket{\text{vac}}$, we can write $\ket{\phi}$ as a
superposition of its $N-$particle components:

$\ket{\phi}=\sum_{M=0}^{\infty}\ket{\phi^{M}}=\sum_{M=0}^{\infty}\sum_{i=1}^{D}\frac{1}{M!}\int_{\omega\in\mathbb{R}^{M}}\phi^{M}_{i}(\omega)\bigg{(}\prod_{j=1}^{M}a_{\omega_{j}}^{\dagger}\bigg{)}\ket{e_{i}}\otimes\ket{\text{vac}},$

where $\\{\ket{e_{i}},i\in\\{1,2\dots D\\}\\}$ is a basis for the
$D-$dimensional local system and

$\phi^{M}_{i}(\omega)=\bra{e_{i}}\otimes\bra{\text{vac}}\prod_{i=1}^{M}a_{\omega_{i}}\ket{\phi}=\bra{e_{i}}\otimes\bra{\text{vac}}U(s_{0},0)\mathcal{T}\bigg{[}\prod_{j=1}^{M}a_{\omega_{j}}(s_{0})\prod_{j=2}^{i}L(s_{j})\bigg{]}\ket{\sigma}\otimes\ket{\text{vac}},$

where $a_{\omega}(s)=U(0,s)a_{\omega}U(s,0)$ and we have used $s_{0}\geq
s_{i}\geq s_{i-1}\dots\geq s_{1}$ to introduce the time-ordering operator. We
then use the Heisenberg equations of motion for $a_{\omega}(s)$ to express

$a_{\omega}(s_{0})=a_{\omega}(0)e^{-i\omega
s}-iv^{*}(\omega)\int_{0}^{s_{0}}L(\tau)e^{-i\omega(s_{0}-\tau)}d\tau.$

We thus obtain

$\displaystyle\phi_{i}^{M}(\omega)$
$\displaystyle=(-i)^{M}\bigg{(}\prod_{j=1}^{M}v^{*}(\omega_{j})\bigg{)}\int_{\tau\in[0,s]^{M}}\bra{e_{i}}\otimes\bra{\text{vac}}U(s,0)\mathcal{T}\bigg{[}\prod_{j=1}^{M}L(\tau_{j})\prod_{j=2}^{i}L(s_{j})\bigg{]}\ket{\sigma}\otimes\ket{\text{vac}}e^{-i\omega\cdot(s-\tau)}d\tau$
$\displaystyle=(-i)^{M}\bigg{(}\prod_{j=1}^{M}v^{*}(\omega_{j})\bigg{)}\int_{\tau\in[0,s]^{M}}\bra{e_{i}}G_{M+i}([\tau,s];s_{0})\ket{\sigma}e^{-i\omega\cdot(s_{0}-\tau)}d\tau.$

This immediately allows us to bound the norm of the $M-$particle component of
$\ket{\phi}$:

$\left\|{\ket{\phi^{M}}}\right\|^{2}\leq\frac{(2\pi)^{M}\left\|{v}\right\|_{\infty}^{2M}}{M!}\int_{\tau\in[0,s]^{M}}\bra{\sigma}G_{M+i}^{\dagger}([\tau,s_{2}\dots
s_{i}];s)G_{M+i}([\tau,s_{2}\dots
s_{i}];s)\ket{\sigma}d\tau\leq\frac{\left\|{v}\right\|_{\infty}^{2M}(2\pi
s)^{M}}{M!}\left\|{L}\right\|_{2}^{2(M+i)}.$

Finally, noting that the restriction of the operator
$\int_{-\infty}^{\infty}v(\omega)a_{\omega}d\omega$ on the $M-$particle
subspace is
$\sqrt{M}\big{[}\int_{-\infty}^{\infty}|v(\omega)|^{2}d\omega\big{]}^{1/2}$,
we obtain that

$\displaystyle\left\|{\bigg{(}\int_{-\infty}^{\infty}v(\omega)a_{\omega}d\omega\bigg{)}U(s,0)L(s_{i-1})\dots
L(s_{2})\ket{\sigma}\otimes\ket{\text{vac}}}\right\|^{2}$
$\displaystyle\leq\bigg{(}\int_{-\infty}^{\infty}|v(\omega)|^{2}d\omega\bigg{)}\sum_{M=0}^{\infty}M\left\|{\ket{\phi^{M}}}\right\|^{2}$
$\displaystyle\leq\left\|{L}\right\|_{2}^{2(i+1)}2\pi
s\left\|{v}\right\|_{\infty}^{2}e^{2\pi\left\|{v}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}s},$

and consequently we obtain that
$\forall\sigma\in\mathcal{H}_{S},s_{0}\in(s_{i},s_{i+1})$,

$\left\|{\partial_{s_{0}}G_{N}(\\{s_{0},s\\};t)\ket{\sigma}}\right\|_{2}\leq\left\|{L}\right\|_{2}^{N}\bigg{(}\sup_{s_{0}\in[0,t]}\frac{\left\|{[L,H_{S}(s_{0})]}\right\|_{2}}{\left\|{L}\right\|_{2}}+\left\|{[L,L^{\dagger}]}\right\|_{2}\sqrt{2\pi
t}\left\|{v}\right\|_{\infty}e^{\pi\left\|{v}\right\|_{\infty}^{2}\left\|{L}\right\|_{2}^{2}t/2}\int_{-\infty}^{\infty}|v(\omega)|^{2}d\omega\bigg{)}.$

Since this bound holds for every $i\in\\{0,1\dots N-2\\}$, it proves that
assumption 2 holds for square integrable coupling functions $\square$.

##  Appendix D Proof of Lemma 7

Proof: Our goal is to develop a Lorentzian approximation to a function
$\Gamma$ which is bounded and differentiable. We develop such an approximation
in two steps:

  1. 1.

We first use the fact that Lorentzians are mollifiers to construct an
approximation $\Gamma_{1}:\mathbb{R}\to\mathbb{R}$ to $\Gamma$ within the
frequency window $[-\omega_{c},\omega_{c}]$

$\displaystyle\Gamma_{1}(\omega)=\frac{\kappa}{\pi}\int_{-\omega_{c}}^{\omega_{c}}\frac{\Gamma(\nu)}{(\omega-\nu)^{2}+\kappa^{2}}d\nu$
(24)

Note that as $\kappa\to 0$, $\Gamma_{1}(\omega)$ is expected to converge to
$\Gamma(\omega)$ for $\omega\in(-\omega_{c},\omega_{c})$. For a non-zero value
of $\kappa$, an error will be incurred in approximating $\Gamma$ by
$\Gamma_{1}$.

  2. 2.

Next, we numerically discretize the integral in Eq. 24 to obtain the
approximation $\hat{\Gamma}:\mathbb{R}\to\mathbb{R}$ to $\Gamma$. 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
$\delta=2\omega_{c}/(M-1)$ and thus

$\displaystyle\hat{\Gamma}(\omega)=\sum_{n=1}^{M}\frac{\kappa\delta}{\pi}\frac{\Gamma(\omega_{n})}{(\omega-\omega_{n})^{2}+\kappa^{2}}.$
(25)

We now consider the error incurred in approximating $\Gamma$ within a
frequency window $[-\omega_{c},\omega_{c}]$ by $\hat{\Gamma}$. The relevant
error of interest is the $\text{L}^{1}$ norm error between
$\Gamma\cdot\text{rect}_{\omega_{c}}$ and $\hat{\Gamma}$, which can be upper
bounded by

$\int_{-\infty}^{\infty}\bigg{|}\Gamma(\omega)\text{rect}_{\omega_{c}}(\omega)-\hat{\Gamma}(\omega)\bigg{|}d\omega\leq\underbrace{\int_{-\omega_{c}}^{\omega_{c}}\big{|}\Gamma(\omega)-\Gamma_{1}(\omega)\big{|}d\omega}_{\textnormal{Error
between $\Gamma$ and
$\Gamma_{1}$}}+\underbrace{\int_{-\omega_{c}}^{\omega_{c}}\big{|}\Gamma_{1}(\omega)-\hat{\Gamma}(\omega)\big{|}d\omega}_{\textnormal{Error
between $\Gamma_{1}$ and
$\hat{\Gamma}$}}+\underbrace{\int_{|\omega|\geq\omega_{c}}|\hat{\Gamma}(\omega)|d\omega}_{\textnormal{Error
due to tails}}.$

We can estimate each of these terms individually.

_Error between $\Gamma$ and $\Gamma_{1}$_: It follows from the definition of
$\Gamma_{1}$ that

$\displaystyle\big{|}\Gamma_{1}(\omega)-\Gamma(\omega)\big{|}$
$\displaystyle=\frac{\kappa}{\pi}\bigg{|}\int_{-\omega_{c}}^{\omega_{c}}\frac{\Gamma(\nu)}{(\omega-\nu)^{2}+\kappa^{2}}d\nu-\int_{-\infty}^{\infty}\frac{\Gamma(\omega)}{(\omega-\nu)^{2}+\kappa^{2}}d\nu\bigg{|}$
$\displaystyle\leq\frac{\gamma}{\pi}\int_{|\nu|\geq\omega_{c}}\bigg{|}\frac{\Gamma(\omega)}{(\omega-\nu)^{2}+\kappa^{2}}\bigg{|}d\nu+\frac{\kappa}{\pi}\int_{-\omega_{c}}^{\omega_{c}}\bigg{|}\frac{\Gamma(\omega)-\Gamma(\nu)}{(\omega-\nu)^{2}+\kappa^{2}}\bigg{|}d\nu.$

The first integral can analytically integrated to obtain:

$\displaystyle\frac{\gamma}{\pi}\int_{|\nu|\geq\omega_{c}}\bigg{|}\frac{\Gamma(\omega)}{(\omega-\nu)^{2}+\kappa^{2}}\bigg{|}d\nu=\Gamma(\omega)\bigg{(}1-\frac{1}{\pi}\tan^{-1}\frac{|\omega-\omega_{c}|}{\kappa}-\frac{1}{\pi}\tan^{-1}\frac{|\omega+\omega_{c}|}{\kappa}\bigg{)}.$

To estimate the second integral, we make use of Taylor’s theorem:
$\forall\nu\in[-\omega_{c},\omega_{c}]:\
|\Gamma(\omega)-\Gamma(\nu)|\leq\Gamma_{\textnormal{max}}^{d}(\omega_{c})|\omega-\nu|$
where
$\Gamma_{\textnormal{max}}^{d}(\omega_{c})=\sup_{\omega\in[-\omega_{c},\omega_{c}]}|\Gamma^{\prime}(\omega)|$.
It then follows that $\forall\omega\in[-\omega_{c},\omega_{c}]$ that

$\displaystyle\frac{\kappa}{\pi}\int_{-\omega_{c}}^{\omega_{c}}\bigg{|}\frac{\Gamma(\omega)-\Gamma(\nu)}{(\omega-\nu)^{2}+\kappa^{2}}\bigg{|}d\nu\leq\frac{\kappa\Gamma_{\text{max}}^{d}(\omega_{c})}{2\pi}\bigg{(}\log\frac{(\omega-\omega_{c})^{2}+\kappa^{2}}{\kappa^{2}}+\log\frac{(\omega+\omega_{c})^{2}+\kappa^{2}}{\kappa^{2}}\bigg{)}.$

We thus obtain that

$\int_{-\omega_{c}}^{\omega_{c}}|\Gamma_{1}(\omega)-\Gamma(\omega)|d\omega\leq\omega_{c}\Gamma_{\text{max}}(\omega_{c})\xi\bigg{(}\frac{\kappa}{\omega_{c}}\bigg{)}+\omega_{c}^{2}\Gamma_{\text{max}}^{d}\xi^{d}\bigg{(}\frac{\kappa}{\omega_{c}}\bigg{)}$

where

$\displaystyle\xi(x)=\frac{1}{\pi}\tan^{-1}\bigg{(}\frac{x}{2}\bigg{)}-\frac{x}{\pi}\log\bigg{(}1+\frac{4}{x^{2}}\bigg{)},$
$\displaystyle\xi^{d}(x)=4x\log\bigg{(}1+\frac{4}{x^{2}}\bigg{)}-4x+4x\tan^{-1}\bigg{(}\frac{2}{x}\bigg{)}.$

We note that $\xi,\xi_{d}=\mathcal{O}(x\text{log}(1/x))$ and therefore

$\int_{-\omega_{c}}^{\omega_{c}}\big{|}\Gamma_{1}(\omega)-\Gamma(\omega)\big{|}d\omega\leq\mathcal{O}\bigg{(}\kappa\Gamma_{\text{max}}\log\frac{\kappa}{\omega_{c}}\bigg{)}+\mathcal{O}\bigg{(}\omega_{c}\kappa\Gamma_{\text{max}}^{d}(\omega_{c})\log\frac{\kappa}{\omega_{c}}\bigg{)}$

_Error between $\Gamma_{1}$ and $\hat{\Gamma}$_: Next, to obtain
$\hat{\Gamma}$ from $\Gamma_{1}$, we simply discretize the integral in the
definition of $\Gamma_{1}$ — the Taylor’s theorem can again be used to bound
the resulting discretization error. In particular,

$\int_{-\omega_{c}}^{\omega_{c}}\big{|}\Gamma_{1}(\omega)-\hat{\Gamma}(\omega)\big{|}d\omega\leq\frac{2\omega_{c}^{3}}{\pi(M-1)}\bigg{(}\frac{\Gamma_{\text{max}}^{d}(\omega_{c})}{\kappa}+\frac{3\sqrt{3}}{8}\frac{\Gamma_{\text{max}}(\omega_{c})}{\kappa^{2}}\bigg{)}=\mathcal{O}\bigg{(}\frac{\omega_{c}^{3}\Gamma_{\text{max}}^{d}(\omega_{c})}{\kappa
M}\bigg{)}+\mathcal{O}\bigg{(}\frac{\omega_{c}^{3}\Gamma_{\text{max}}}{\kappa^{2}M}\bigg{)}.$

_Error due to tails:_ This error can be explicitly evaluated.

$\displaystyle\int_{|\omega|\geq\omega_{c}}|\hat{\Gamma}(\omega)|d\omega$
$\displaystyle\leq\sum_{n=1}^{M}\frac{\kappa\delta\Gamma_{\text{max}}}{\pi}\int_{|\omega|\geq\omega_{c}}\frac{1}{(\omega-\omega_{n})^{2}+\kappa^{2}}d\omega,$
$\displaystyle=\sum_{n=1}^{M}\frac{\delta\Gamma_{\text{max}}}{\pi}\bigg{(}\tan^{-1}\bigg{(}\frac{\kappa}{\omega_{c}-\omega_{n}}\bigg{)}+\tan^{-1}\bigg{(}\frac{\kappa}{\omega_{c}+\omega_{n}}\bigg{)}\bigg{)},$
$\displaystyle=\delta\Gamma_{\text{max}}+\frac{2\delta\Gamma_{\text{max}}}{\pi}\sum_{n=1}^{M-1}\tan^{-1}\bigg{(}\frac{\kappa}{n\delta}\bigg{)}.$

Noting that for $x>0$, $\tan^{-1}(1/x)$ is convex upwards, it follows that for
any $\theta>0$

$\sum_{n=1}^{M-1}\tan^{-1}\bigg{(}\frac{1}{n\theta}\bigg{)}\leq\frac{1}{\theta}\int_{0}^{(M-1)\theta}\tan^{-1}\bigg{(}\frac{1}{x}\bigg{)}dx=(M-1)\tan^{-1}\bigg{(}\frac{1}{(M-1)\theta}\bigg{)}+\frac{1}{2\theta}\log\big{(}1+(M-1)^{2}\theta^{2}\big{)}.$

Consequently,

$\int_{|\omega|\geq\omega_{c}}|\hat{\Gamma}(\omega)|d\omega\leq\mathcal{O}\bigg{(}\frac{\Gamma_{\text{max}}}{M}\bigg{)}+\mathcal{O}\bigg{(}\kappa\Gamma_{\text{max}}\log\frac{\kappa}{\omega_{c}}\bigg{)}.$

This proves the lemma $\square$.

##  Appendix E Proof of lemma 9

Proof: We start by noting that
$H_{\omega_{c}}^{E}=\int_{-\omega_{c}}^{\omega_{c}}\omega
a_{\omega}^{\dagger}a_{\omega}d\omega$ and therefore

$e^{isH_{\omega_{c}}^{E}}Be^{-isH_{\omega_{c}}^{E}}=\int_{-\omega_{c}}^{\omega_{c}}v(\omega)a_{\omega}e^{-i\omega
s}d\omega.$

Evaluating $e^{isH_{\omega_{c},N}^{E}}Be^{-isH_{\omega_{c},N}^{E}}$ is more
involved since we need to diagonalize $H_{\omega_{c},N}^{E}$. It is convenient
to introduce the polynomial $\pi_{i}$, which is $p_{i}$ after normalization

$\pi_{i}(\omega)=\frac{p_{i}(\omega)\big{(}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}d\omega\big{)}^{1/2}}{\big{(}\int_{-\omega_{c}}^{\omega_{c}}p_{i}^{2}(\omega)|v(\omega)|^{2}d\omega\big{)}^{1/2}}$

Suppose $\Phi=\sum_{i=0}^{N-1}\phi_{i}A_{i}$ is an eigenmode at eigenfrequency
$\Omega$ i.e. $[\Phi,H_{\omega_{c},N}^{E}]=\Omega\Phi$. Using Eq. 19, we
obtain

$(\alpha_{i}-\Omega)\phi_{i}+\sqrt{\beta_{i}}\phi_{i-1}+\sqrt{\beta_{i+1}}\phi_{i+1}=0.$

It then follows that for $i\in\\{0,1,2\dots N-1\\}$,
$\phi_{i}=\pi_{i}(\omega)\phi_{0}$ and $p_{N}(\Omega)=0$. Thus the possible
eigenenergies of $H^{E}_{\omega_{c},N}$, denoted by
$\Omega_{0},\Omega_{1}\dots\Omega_{N-1}$, are roots of the $N^{\text{th}}$
orthogonal polynomial. For $i\in\\{0,1,2\dots N-1\\}$, the eigenmode
$\Phi_{i}$ associated with $\Omega_{i}$ can be written as
$\Phi_{i}=\sum_{j=0}^{N-1}\phi^{i}_{j}A_{j}$ where

$\displaystyle\phi_{i}^{j}=\frac{\pi_{i}(\Omega_{j})}{N_{j}}\ \text{where
}N_{j}=\bigg{[}\sum_{i=0}^{N-1}\pi_{i}^{2}(\Omega_{j})\bigg{]}^{1/2}.$ (26)

We note that from the orthogonality of the eigenmodes $\phi_{i}$, it follows
that

$\displaystyle\sum_{i=0}^{N-1}\pi_{i}(\Omega_{j})\pi_{i}(\Omega_{j^{\prime}})=N_{j}^{2}\delta_{j,j^{\prime}}\text{
and
}\sum_{j=0}^{N-1}\pi_{i}(\Omega_{j})\pi_{i^{\prime}}(\Omega_{j})=N_{j}^{2}\delta_{i,i^{\prime}}.$
(27)

Finally, consider evaluating
$e^{isH_{\omega_{c},N}^{E}}Be^{-isH_{\omega_{c},N}^{E}}=\big{(}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}d\omega\big{)}^{1/2}\big{(}e^{isH_{\omega_{c},N}^{E}}A_{0}e^{-isH_{\omega_{c},N}^{E}}\big{)}$.
Noting that $A_{0}=\sum_{i=0}^{N-1}\phi_{0}^{j}\Phi_{j}$, we obtain

$e^{isH_{\omega_{c},N}^{E}}Be^{-isH_{\omega_{c},N}^{E}}=\bigg{(}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}d\omega\bigg{)}^{1/2}\sum_{j=0}^{N-1}\phi_{0}^{j}\Phi_{j}e^{-i\Omega_{j}s}.$

Finally, we evaluate the commutator $[\delta B_{N}(s),\delta
B_{N}^{\dagger}(s)]$ — with straightforward computation, it follows that

$\displaystyle[\delta B_{N}(s),\delta
B_{N}^{\dagger}(s)]=2\int_{-\infty}^{\infty}|v(\omega)|^{2}d\omega-2\sum_{i,j=0}^{N-1}\phi^{j}_{0}\phi^{j}_{i}{\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}\pi_{i}(\omega)\cos((\omega-\Omega_{j})s)d\omega}$
(28)

Next we use the Gauss quadrature theorem, from which it would follow that for
a given $N>0,\ \exists w_{0},w_{1}\dots w_{N-1}>0$ with
$\sum_{i=0}^{N-1}w_{i}=1$ such that for all polynomials $P$ of degree $\leq
2N-1$,

$\int_{-\omega_{c}}^{\omega_{c}}P(\omega)|v(\omega)|^{2}d\omega=\bigg{(}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}d\omega\bigg{)}\sum_{i=0}^{N-1}w_{i}P(\Omega_{i}).$

To use this result, we note that from Taylor’s theorem it follows that the
cosine in Eq. 28 can be expanded into the sum of a degree $N$ polynomial and a
remainder i.e. $\forall j\in\\{0,1,2\dots
N-1\\},\forall\omega\in[-\omega_{c},\omega_{c}]:\cos((\omega-\Omega_{j})s)=C_{N}(\omega-\Omega_{j})+R_{N}(\omega-\Omega_{j})$
where $C_{N}$ is a degree $N$ polynomial with $C_{N}(0)=1$ and $R_{N}$ is the
remainder which satisfies the estimate
$\left\|{R_{N}}\right\|_{\infty}\leq(2\omega_{c}s)^{N+1}/(N+1)!$. With this
decomposition of the cosine function, Eq. 28 reduces to

$\displaystyle[\delta B_{N}(s),\delta B_{N}^{\dagger}(s)]$
$\displaystyle=2\int_{-\infty}^{\infty}|v(\omega)|^{2}d\omega-2\sum_{i,j=0}^{N-1}\phi_{0}^{j}\phi_{i}^{j}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}\pi_{i}(\omega)C_{N}(\omega-\Omega_{j})d\omega-2\sum_{i,j=0}^{N-1}\phi_{0}^{j}\phi_{i}^{j}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}R_{N}(\omega-\Omega_{j})\pi_{i}(\omega)d\omega.$

Noting that for $i\leq N-1$, $\pi_{i}(\omega)C_{N}(\omega-\Omega_{j})$ is a
polynomial of degree $\leq 2N-1$, we obtain by an application of the Gauss
quadrature theorem and using Eq. 26 that

$\displaystyle\sum_{i,j=0}^{N-1}\phi_{0}^{j}\phi_{i}^{j}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}\pi_{i}(\omega)C_{N}(\omega-\Omega_{j})d\omega$
$\displaystyle\qquad=\bigg{(}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}d\omega\bigg{)}\sum_{i,j,k=0}^{N-1}w_{k}\frac{\pi_{i}(\Omega_{j})\pi_{i}(\Omega_{k})}{N_{j}^{2}}C_{N}(\Omega_{k}-\Omega_{j})$
$\displaystyle\qquad=\bigg{(}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}d\omega\bigg{)}\sum_{i,k=0}^{N-1}w_{k}\delta_{i,k}C_{N}(\Omega_{k}-\Omega_{j})=\bigg{(}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}d\omega\bigg{)}\sum_{k=0}^{N-1}w_{k}=\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}d\omega.$

Consequently, it follows that

$[\delta B_{N}(s),\delta
B_{N}^{\dagger}(s)]=-2\sum_{i,j=0}^{N-1}\phi_{0}^{j}\phi_{i}^{j}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}R_{N}(\omega-\Omega_{j})\pi_{i}(\omega)d\omega.$

Noting $|\phi_{i}^{j}|\leq 1$ for $i,j\in\\{0,1,2\dots N-1\\}$, we then obtain
the following estimate

$\displaystyle\Big{|}[\delta B_{N}(s),\delta B_{N}^{\dagger}(s)]\Big{|}$
$\displaystyle\leq
2\sum_{i,j=0}^{N-1}\bigg{|}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}R_{N}(\omega-\Omega_{j})\pi_{i}(\omega)d\omega\bigg{|}$
$\displaystyle\leq
2N\left\|{v}\right\|_{\infty}\left\|{R_{N}}\right\|_{\infty}\sum_{i=0}^{N-1}\bigg{|}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|\pi_{i}(\omega)d\omega\bigg{|}$
$\displaystyle\leq
2N\left\|{v}\right\|_{\infty}\frac{(2\omega_{c}s)^{N+1}}{(N+1)!}\sum_{i=0}^{N-1}\bigg{(}2\omega_{c}\int_{-\omega_{c}}^{\omega_{c}}|v(\omega)|^{2}\pi_{i}^{2}(\omega)d\omega\bigg{)}^{1/2}$
$\displaystyle\leq
4\omega_{c}\left\|{v}\right\|_{\infty}^{2}\frac{(2\omega_{c}s)^{N+1}}{(N-1)!}=\mathcal{O}\bigg{(}\frac{\left\|{v}\right\|_{\infty}^{2}(2\omega_{c})^{N+2}s^{N+1}}{(N-1)!}\bigg{)},$

wherein in the last step we have used the fact that the polynomial
$\pi_{i}(\omega)$ is normalized. This completes the proof $\square$.

## References

  * Breuer  _et al._ (2002) H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Finsterhölzl  _et al._ (2020) R. Finsterhölzl, M. Katzer, and A. Carmele, Physical Review B 102, 174309 (2020).
  * Chin  _et al._ (2011) A. W. Chin, J. Prior, S. F. Huelga, and M. B. Plenio, Physical review letters 107, 160601 (2011).
  * Groeblacher  _et al._ (2015) S. Groeblacher, A. Trubarov, N. Prigge, G. Cole, M. Aspelmeyer, and J. Eisert, Nature communications 6, 1 (2015).
  * de Vega  _et al._ (2008) I. de Vega, D. Porras, and J. I. Cirac, Physical review letters 101, 260404 (2008).
  * Calajó  _et al._ (2019) G. Calajó, Y.-L. L. Fang, H. U. Baranger, F. Ciccarello, _et al._ , Physical review letters 122, 073601 (2019).
  * Andersson  _et al._ (2019) G. Andersson, B. Suri, L. Guo, T. Aref, and P. Delsing, Nature Physics 15, 1123 (2019).
  * González-Tudela  _et al._ (2019) A. González-Tudela, C. S. Muñoz, and J. I. Cirac, Physical review letters 122, 203603 (2019).
  * Aref  _et al._ (2016) T. Aref, P. Delsing, M. K. Ekström, A. F. Kockum, M. V. Gustafsson, G. Johansson, P. J. Leek, E. Magnusson, and R. Manenti, in  _Superconducting devices in quantum optics_ (Springer, 2016) pp. 217–244.
  * Leonforte  _et al._ (2021) L. Leonforte, A. Carollo, and F. Ciccarello, Physical Review Letters 126, 063601 (2021).
  * Ishizaki and Tanimura (2005) A. Ishizaki and Y. Tanimura, Journal of the Physical Society of Japan 74, 3131 (2005).
  * Chin  _et al._ (2012) A. Chin, S. F. Huelga, and M. B. Plenio, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 3638 (2012).
  * Ivanov and Breuer (2015) A. Ivanov and H.-P. Breuer, Physical Review A 92, 032113 (2015).
  * Caycedo-Soler  _et al._ (2021) F. Caycedo-Soler, A. Mattioni, J. Lim, T. Renger, S. Huelga, and M. Plenio, arXiv preprint arXiv:2106.14286 (2021).
  * Xu  _et al._ (2018) M. Xu, Y. Yan, Y. Liu, and Q. Shi, The Journal of chemical physics 148, 164101 (2018).
  * Smirne and Vacchini (2010) A. Smirne and B. Vacchini, Physical Review A 82, 022110 (2010).
  * Pereverzev and Bittner (2006) A. Pereverzev and E. R. Bittner, The Journal of chemical physics 125, 104906 (2006).
  * Kidon  _et al._ (2015) L. Kidon, E. Y. Wilner, and E. Rabani, The Journal of chemical physics 143, 234110 (2015).
  * Vacchini and Breuer (2010) B. Vacchini and H.-P. Breuer, Physical Review A 81, 042103 (2010).
  * Schröder  _et al._ (2007) M. Schröder, M. Schreiber, and U. Kleinekathöfer, The Journal of chemical physics 126, 114102 (2007).
  * Timm (2011) C. Timm, Physical Review B 83, 115416 (2011).
  * Mukamel  _et al._ (1978) S. Mukamel, I. Oppenheim, and J. Ross, Physical Review A 17, 1988 (1978).
  * Albeverio  _et al._ (2007) S. Albeverio, L. Cattaneo, S. Mazzucchi, and L. Di Persio, Journal of mathematical physics 48, 102109 (2007).
  * Makarov and Makri (1994) D. E. Makarov and N. Makri, Chemical physics letters 221, 482 (1994).
  * Smith and Caldeira (1987) C. M. Smith and A. Caldeira, Physical Review A 36, 3509 (1987).
  * Whalen  _et al._ (2017) S. Whalen, A. Grimsmo, and H. Carmichael, Quantum Science and Technology 2, 044008 (2017).
  * Grimsmo (2015) A. L. Grimsmo, Physical review letters 115, 060402 (2015).
  * Jørgensen and Pollock (2019) M. R. Jørgensen and F. A. Pollock, Physical review letters 123, 240602 (2019).
  * Del Pino  _et al._ (2018) J. Del Pino, F. A. Schröder, A. W. Chin, J. Feist, and F. J. Garcia-Vidal, Physical review letters 121, 227401 (2018).
  * Strathearn  _et al._ (2018) A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett, Nature communications 9, 1 (2018).
  * Rosenbach  _et al._ (2016) R. Rosenbach, J. Cerrillo, S. F. Huelga, J. Cao, and M. B. Plenio, New Journal of Physics 18, 023035 (2016).
  * Huh  _et al._ (2014) J. Huh, S. Mostame, T. Fujita, M.-H. Yung, and A. Aspuru-Guzik, New Journal of Physics 16, 123008 (2014).
  * Chin  _et al._ (2010) A. W. Chin, Á. Rivas, S. F. Huelga, and M. B. Plenio, Journal of Mathematical Physics 51, 092109 (2010).
  * Woods  _et al._ (2014) M. Woods, R. Groux, A. Chin, S. F. Huelga, and M. B. Plenio, Journal of Mathematical Physics 55, 032101 (2014).
  * Fang (2019) Y.-L. L. Fang, Computer Physics Communications 235, 422 (2019).
  * Vidal (2004) G. Vidal, Physical review letters 93, 040502 (2004).
  * Daley  _et al._ (2004) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Journal of Statistical Mechanics: Theory and Experiment 2004, P04005 (2004).
  * Banuls  _et al._ (2009) M.-C. Banuls, M. B. Hastings, F. Verstraete, and J. I. Cirac, Physical review letters 102, 240603 (2009).
  * Daley  _et al._ (2009) A. Daley, J. Taylor, S. Diehl, M. Baranov, and P. Zoller, Physical review letters 102, 040402 (2009).
  * Verstraete  _et al._ (2004) F. Verstraete, J. J. Garcia-Ripoll, and J. I. Cirac, Physical review letters 93, 207204 (2004).
  * Berry  _et al._ (2014) D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, in  _Proceedings of the forty-sixth annual ACM symposium on Theory of computing_ (2014) pp. 283–292.
  * Low and Chuang (2019) G. H. Low and I. L. Chuang, Quantum 3, 163 (2019).
  * Berry  _et al._ (2015) D. W. Berry, A. M. Childs, and R. Kothari, in  _2015 IEEE 56th Annual Symposium on Foundations of Computer Science_ (IEEE, 2015) pp. 792–809.
  * Zanardi  _et al._ (2016) P. Zanardi, J. Marshall, and L. C. Venuti, Physical Review A 93, 022312 (2016).
  * Chenu  _et al._ (2017) A. Chenu, M. Beau, J. Cao, and A. del Campo, Physical review letters 118, 140403 (2017).
  * Cleve and Wang (2016) R. Cleve and C. Wang, arXiv preprint arXiv:1612.09512 (2016).
  * Budini (2000) A. A. Budini, Physical Review A 63, 012106 (2000).
  * Megier  _et al._ (2018) N. Megier, W. T. Strunz, C. Viviescas, and K. Luoma, Physical review letters 120, 150402 (2018).
  * Tamascelli  _et al._ (2018) D. Tamascelli, A. Smirne, S. F. Huelga, and M. B. Plenio, Physical review letters 120, 030402 (2018).
  * Pleasance  _et al._ (2020) G. Pleasance, B. M. Garraway, and F. Petruccione, Physical Review Research 2, 043058 (2020).
  * Mazzola  _et al._ (2009) L. Mazzola, S. Maniscalco, J. Piilo, K.-A. Suominen, and B. M. Garraway, Physical Review A 80, 012104 (2009).
  * Dalton  _et al._ (2001) B. Dalton, S. M. Barnett, and B. Garraway, Physical Review A 64, 053813 (2001).
  * Garraway and Dalton (2006) B. Garraway and B. Dalton, Journal of Physics B: Atomic, Molecular and Optical Physics 39, S767 (2006).
  * Dalton  _et al._ (2012) B. Dalton, S. M. Barnett, and B. Garraway, in  _Coherence and Quantum Optics VIII: Proceedings of the Eighth Rochester Conference on Coherence and Quantum Optics, held at the University of Rochester, June 13–16, 2001_ (Springer Science & Business Media, 2012) p. 495.
  * Fang  _et al._ (2018) Y.-L. L. Fang, F. Ciccarello, and H. U. Baranger, New Journal of Physics 20, 043035 (2018).
  * Fan  _et al._ (2010) S. Fan, Ş. E. Kocabaş, and J.-T. Shen, Physical Review A 82, 063821 (2010).
  * Barut and Kraus (1983) A. Barut and J. Kraus, Foundations of Physics 13, 189 (1983).
  * Garraway (1997) B. Garraway, Physical Review A 55, 2290 (1997).
  * Cilluffo  _et al._ (2020) D. Cilluffo, A. Carollo, S. Lorenzo, J. A. Gross, G. M. Palma, and F. Ciccarello, Physical Review Research 2, 043070 (2020).
  * Van Dooren (1995) P. Van Dooren, in  _Short Course, Benelux Meeting on Systems and Control_ , Vol. 120 (Citeseer, 1995).
  * de Vega  _et al._ (2015) I. de Vega, U. Schollwöck, and F. A. Wolf, Physical Review B 92, 155126 (2015).